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Abstract 


Large-eddy  simulation  (LES)  requires  very  high  resolution  in  high  Reynolds  num¬ 
ber,  attached  turbulent  boundary  layers  due  to  the  need  to  capture  the  small, 
dynamically  important  near-wall  eddies.  Resolving  these  eddies  causes  the  compu¬ 
tational  expense  of  the  LES  to  scale  almost  as  strongly  with  the  Reynolds  number 
as  direct  numerical  simulation  for  these  flows.  Wall  modeling  is  a  technique  which 
enables  LES  to  be  performed  on  grids  that  do  not  resolve  the  wall  layer.  Instead, 
it  provides  approximate  boundary  conditions  to  the  LES  at  solid  boundaries,  thus 
allowing  a  much  weaker  scaling  of  the  LES  grid  size  with  the  Reynolds  number. 

Unfortunately,  wall  models  based  on  purely  physical  reasoning  often  lead  to  an 
inaccurate  LES,  particularly  on  coarse  grids  and  at  high  Reynolds  numbers,  be¬ 
cause  they  do  not  account  for  the  numerical  and  SGS  modeling  errors  that  become 
large  in  these  types  of  simulations.  To  address  these  errors,  optimal  control-based 
wall  models  have  been  developed  by  previous  investigators.  While  these  have  the 
demonstrated  ability  to  account  for  the  aforementioned  errors,  they  have  two  pri¬ 
mary  drawbacks:  1)  high  computational  expense,  due  to  the  optimization  proce¬ 
dure,  and  2)  a  lack  of  predictability,  because  the  control  targets  are  prescribed  a 
priori. 

The  goal  of  this  work  is  to  address  these  two  issues  in  order  to  make  control- 
based  wall  modeling  feasible  for  engineering  applications.  To  reduce  the  expense, 
the  adjoint  equations,  which  are  used  to  determine  the  gradients  needed  for  the 
optimization,  have  been  reformulated  to  minimize  the  effort  required  in  the  op¬ 
timization  procedure.  Further,  the  optimization  algorithm  has  been  modified  to 
only  use  near-wall  information  so  no  work  is  wasted  in  regions  of  the  flow  which 

iii 


are  insensitive  to  the  control.  Such  an  approach  reduces  the  computational  cost  of 
the  method  by  an  order  of  magnitude  without  a  reduction  in  the  accuracy  of  the 
simulation. 

To  make  the  method  predictive,  a  near-wall  Reynolds-averaged  Navier-Stokes 
(RANS)  model  has  been  coupled  to  the  LES/controller  system  to  provide  a  target 
for  the  control.  This  coupling  is  accomplished  by  using  the  LES  to  provide  the 
velocity  boundary  conditions  for  RANS  away  from  the  wall,  while  the  RANS  feeds 
back  into  the  LES  through  the  definition  of  the  cost  function  that  is  minimized 
by  the  control.  An  additional  degree  of  coupling  enables  the  RANS  to  provide  the 
mean  wall  stress  for  the  LES.  The  control  then  provides  the  fluctuating  wall  stress 
which  minimizes  the  cost  function.  Using  this  method  in  plane  channel  flow,  an 
accurate  prediction  of  the  mean  velocity  profile  has  been  obtained  over  a  range  of 
Reynolds  numbers  and  on  different  grids.  The  results  are  comparable  to  those  from 
previous  control-based,  non-predictive  models,  and  are  much  more  accurate  than 
the  predictions  of  traditional  wall  models. 
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Chapter  1 
Introduction 


1.1  Wall  Modeling  Background 

The  ability  to  accurately  simulate  fluid  flows  has  important  applications  in  en¬ 
gineering  design  and  analysis.  One  of  the  most  significant  impediments  to  such 
simulations  is  the  change  from  laminar  flow,  in  which  flow  features  are  present  only 
over  a  small  number  of  spatial  and  temporal  scales,  to  turbulent  flow,  where  a  very 
wide  range  of  dynamically  important  scales  in  both  space  and  time  are  present.  Re¬ 
solving  all  these  scales  is  the  most  serious  impediment  to  high  fidelity  simulations 
of  fluid  dynamics.  It  has  been  estimated  that  the  required  number  of  grid  points  for 
a  fully  resolved  simulation  scales  as  Re9/4,  where  Re  is  the  Reynolds  number  which 
measures  the  relative  importance  of  inertial  and  viscous  forces.  A  simulation  that 
resolves  all  flow  scales,  and  hence  requires  no  models,  is  called  a  direct  numerical 
simulation  (DNS).  For  a  recent  review  of  DNS,  see  Moin  and  Mahesh  (1998). 

In  an  effort  to  mitigate  the  high  computational  expense  associated  with  DNS, 
the  technique  of  large-eddy  simulation  (LES)  has  been  developed.  The  computa¬ 
tional  cost  is  reduced  by  applying  a  low-pass  filter  to  the  turbulent  flow,  thereby 
eliminating  many  of  the  small  scales  from  the  LES  field.  From  a  physical  and  engi¬ 
neering  perspective,  the  high  frequency  information  that  is  lost  tends  to  be  of  less 
importance  to  practical  problems.  However,  the  short  wavelength  physics  can  have 
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CHAPTER  1.  INTRODUCTION 


a  significant  impact  on  the  evolution  of  the  flow,  and  so  its  effects  on  the  LES  field 
are  incorporated  through  the  use  of  models.  These  models  are  denoted  sub-grid 
scale  (SGS)  models  to  indicate  that  they  supply  information  from  scales  too  small 
to  be  accurately  captured  by  the  numerical  grid.  Much  effort  has  been  put  into 
developing  effective  models  and  techniques  to  perform  LES,  and  good  introductions 
to  these  and  other  issues  found  in  LES  are  provided  by  Carati  (2001),  Sagaut  (2002) 
and  Meneveau  and  Katz  (2000). 

Over  the  years,  models  have  been  developed  that  allow  LES  to  be  successfully 
applied  in  many  types  of  flow  situations.  One  area,  however,  that  has  proved 
particularly  challenging  for  SGS  models  is  the  near-wall  region  of  attached  flows. 
This  is  primarily  due  to  the  fact  that  near  the  wall,  flow  structures  scale  in  viscous 
units.  Hence,  if  the  grid  spacing  is  set  to  capture  the  large-,  or  integral-length, 
scales  of  the  flow,  then  near  the  wall,  many  of  the  important  physical  scales  of  the 
flow  become  small  relative  to  the  grid.  In  addition,  flow  structures  in  this  area  tend 
to  be  anisotropic,  and  since  SGS  models  are  designed  to  model  isotropic  eddies  that 
represent  only  a  small  fraction  of  the  total  energy  of  the  flow,  they  cannot  accurately 
represent  the  turbulent  stresses  in  the  vicinity  of  a  wall  (Jimenez  and  Moser,  2000). 
The  number  of  grid  points  required  to  resolve  the  near  wall  shear  stress  producing 
eddies  scales  as  Re2r  (Baggett  et  al.,  1997).  This  makes  the  near-wall  resolution 
requirements  of  LES  almost  as  high  as  DNS. 

In  order  to  perform  simulations  of  attached  flows  at  high  Reynolds  numbers, 
wall  models  have  been  introduced  to  supply  boundary  conditions  to  the  LES  in 
an  effort  to  eliminate  the  need  to  resolve  the  features  near  the  wall.  This  is  the 
reason  the  use  of  wall  modeling  in  LES  is  almost  as  old  as  LES  itself  (Deardorff, 
1970;  Schumann,  1975):  the  computational  expense  when  the  near-wall  region  is 
not  resolved  becomes  much  more  tractable.  Wall  models  continue  to  be  of  interest 
to  this  day  because  of  the  desire  to  simulate  flows  at  the  high  Reynolds  numbers 
found  in  many  engineering  applications.  Examples  of  many  efforts  in  this  field  can 
be  found  in  the  reviews  of  Piomelli  and  Balaras  (2002)  and  Cabot  and  Moin  (2000). 
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A  typical  wall  model  is  one  that  replaces  the  standard  no-slip  velocity  bound¬ 
ary  conditions  at  a  solid  surface  with  approximate  conditions  to  enable  the  LES  to 
accurately  capture  the  large-scale  features  of  the  flow  away  from  the  surface  with¬ 
out  the  inner  layer  being  resolved.  In  addition  to  the  strong  wall-normal  velocity 
gradients,  this  region  also  contains  many  streaky  structures  that  scale  in  the  inner 
units.  The  structures  are  known  to  be  important  for  the  generation  and  transport 
of  turbulent  kinetic  energy  and  shear  stress.  A  fully  resolved  LES  must  resolve  the 
bulk  of  these  features.  Approximate  boundary  conditions  instead  account  for  the 
effects  of  the  near-wall  turbulence  on  the  outer  flow. 

One  set  of  approximate  boundary  conditions  that  have  several  advantages  are 
wall  stresses.  They  are  directly  related  to  the  large  scale  body  forces  and  acceler¬ 
ations  present  in  the  flow  since  they  are  some  of  the  few  external  forces  that  can 
act  on  the  fluid.  This  relationship  implies  that  they  must  be  correct  for  the  flow 
to  be  accurate.  In  addition,  it  is  possible  to  relate  the  stresses  directly  to  the  state 
of  the  flow  in  the  mean  sense  through  a  known  mean  velocity  profile.  As  will  be 
shown  in  subsequent  sections  in  this  chapter,  many  methods  have  been  developed 
that  utilize  such  a  relationship. 

1.1.1  Alternatives  to  Wall  Models 

OfF-the-Wall  Models 

Before  giving  a  detailed  overview  of  wall  models,  it  will  be  useful  to  consider  alter¬ 
natives  to  them  to  motivate  their  necessity.  The  first  alternative  to  be  considered 
is  the  use  of  off  the  wall  Dirichlet  boundary  conditions.  This  type  of  method  cuts 
the  LES  off  above  the  wall  layer  so  there  is  no  need  to  simulate  the  near-wall  re¬ 
gion.  Instead,  velocities  are  prescribed  where  the  LES  is  cut-off,  and  the  simulation 
is  performed  normally  otherwise.  Using  these  boundary  conditions,  it  is  possible 
to  utilize  a  grid  designed  to  capture  the  outer  scales  of  the  flow.  Unfortunately, 
prescription  of  these  velocities  can  prove  challenging,  as  demonstrated  by  Baggett 
et  al.  (1997). 

In  Baggett  et  al.,  a  resolved  LES  was  performed  and  the  velocity  history  at  the 
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cut-off  region  was  recorded.  This  velocity  history  was  then  used  directly  as  an  off 
the  wall  boundary  condition  for  an  LES  that  did  not  resolve  the  wall,  successfully 
recovering  the  resolved  LES  solution.  While  this  demonstrated  the  theoretical 
feasibility  of  this  approach,  difficulties  were  encountered  when  more  challenging 
tests  were  attempted.  Next,  the  velocity  history  was  distorted  while  maintaining  a 
constant  energy  level  to  test  the  sensitivity  of  the  simulation  to  the  the  boundary 
conditions.  When  the  phase  of  the  boundary  data  was  scrambled  but  retained  the 
same  spectra  and  cospectra,  the  simulation  was  still  able  to  produce  reasonably 
accurate  results.  However,  higher  levels  of  scrambling  that  disrupted  these  spectra, 
and  only  retained  the  second-order  statistics,  created  an  artificial  buffer  layer  above 
the  cut-off  layer  before  the  flow  transitioned  to  a  logarithmic  profile.  This  result 
demonstrated  the  need  for  a  significant  amount  of  physical  information,  including 
turbulent  fluctuations,  to  be  included  in  any  off-wall  boundary  conditions. 

The  approach  of  Baggett  was  extended  by  Nicoud  et  al.  (1998)  and  Jimenez  and 
Vasco  (1998)  with  similar  results.  The  former  group  used  a  scaled  velocity  from  the 
interior  of  the  flow  as  the  boundary  condition.  This  was  done  by  assuming  that  the 
velocities  at  two  wall-parallel  planes  had  self  similar  time  scales  so  that  the  two  could 
be  related.  It  was  determined  the  scaling  ratio  needed  to  be  determined  dynamically 
from  the  LES  to  obtain  the  best  results.  With  this  done,  the  statistics  remained 
symmetric  across  the  channel  despite  the  fact  that  this  boundary  condition  was  only 
applied  on  one  side  of  the  channel  while  a  no-slip  boundary  condition  was  used  at 
the  other.  Unfortunately,  when  this  boundary  condition  was  used  at  both  walls,  it 
was  found  that  the  accuracy  of  the  LES  was  diminshed.  No  definitive  conclusions 
could  be  drawn,  however  the  authors  suggested  that  having  one  physical  boundary 
condition  helped  to  maintain  a  physical  realization  of  the  flow. 

The  work  of  Jimenez  and  Vasco  (1998)  involved  prescribing  velocity  boundary 
conditions  at  the  center  of  the  channel  as  a  feasibility  study  before  attempting  to 
prescribe  velocities  in  the  energetically  and  dynamically  important  wall  layer.  When 
scrambled  velocity  data  from  a  full  channel  were  provided  as  boundary  conditions, 
results  similar  to  those  of  Baggett  (1997)  were  obtained.  In  an  effort  to  design  a 
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predictive  model,  the  velocities  from  a  plane  near  the  boundary  at  the  previous 
time  step  were  used  as  a  boundary  condition,  after  being  scaled  to  match  the 
known  rms  fluctuations  at  the  center.  The  velocities  were  further  modified  to 
ensure  uv  =  vw  =  0.  However,  this  produced  an  unphysical  peak  in  the  pressure 
fluctuations  near  the  upper  boundary.  This  phenomenon  was  somewhat  mitigated 
by  setting  the  transpiration  velocity  to  satisfy  continuity  requirements  based  on  the 
gradients  of  the  other  two  velocity  components,  but  in  the  end  the  results  were  not 
accurate  enough  to  warrant  further  investigation. 

Another  approach  to  off-the-wall  boundary  conditions  was  recently  proposed 
by  Iovieno  et  al.  (2004).  In  this  work,  it  was  noted  that  if  the  filter  size  does  not 
decrease  to  zero  as  the  distance  to  the  wall  becomes  small,  then  the  unfiltered  no-slip 
conditions  no  longer  apply.  The  filter  width  is  then  taken  as  a  function  of  the  wall- 
normal  distance  with  a  minimum  size  such  that  the  near  wall  structures  can  still 
be  resolved.  By  expanding  the  velocity  near  the  wall  in  a  Taylor  series,  and  doing 
likewise  for  the  filter  width,  the  corresponding  boundary  conditions  off  the  wall  can 
be  obtained.  However,  due  to  the  need  for  an  accurate  expansion  of  the  variables, 
the  off  the  wall  boundary  condition  must  be  imposed  between  1  <  y+  <  7.  While 
the  method  produces  reasonable  results  at  low  to  moderate  Reynolds  numbers,  the 
proximity  of  the  boundary  conditions  to  the  wall  require  that  most  of  the  near- wall 
turbulence  be  resolved.  Hence,  this  method  is  best  viewed  as  a  means  of  correcting 
boundary  conditions  in  wall-resolved  LES  for  a  non-zero  filter  width  at  the  wall 
rather  than  a  technique  that  will  extend  LES  to  very  high  Reynolds  numbers. 


Hybrid  RANS/LES  Approaches 

A  second  alternative  to  wall  modeling  involves  merging  LES  and  RANS  directly 
into  a  hybrid  simulation.  Since  LES  requires  high  resolution  near  the  wall,  RANS 
equations  are  instead  used  in  this  region  to  reduce  the  number  of  grid  points.  This 
is  because  RANS  eddy  viscosity  models  are  designed  to  supply  all  of  the  turbulent 
stress,  as  opposed  to  LES  SGS  models  which  provide  only  a  small  fraction  thereof. 
Thus,  a  RANS  layer  is  used  as  part  of  the  simulation  near  the  wall  to  account  for 
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more  unresolved  stress.  This  means  the  no-slip  boundary  conditions  can  be  directly 
applied. 

The  difficulty  with  the  hybrid  approach  comes  from  providing  the  matching 
conditions  at  the  boundary  between  the  two  simulations.  The  LES  requires  a 
fluctuating  field  that  transports  turbulent  stresses  across  the  interface,  while  the 
RANS  can  only  provide  a  mean  field  without  the  turbulent  fluctuations.  Such  a 
model  was  suggested  by  Quemere  et  al.  (2000),  who  attempted  to  resolve  this  issue 
by  using  either  a  predictor  simulation  or  by  adding  random  perturbations  to  the 
mean  field.  In  these  cases,  the  method  and  results  become  quite  similar  to  the 
off-the-wall  boundary  conditions  used  by  Baggett  (1997)  and  Nicoud  et  al.  (1998). 

Alternatively,  SGS  models  have  been  developed  that  behave  like  RANS  models 
near  the  wall,  allowing  this  region  to  be  resolved  only  to  the  degree  required  for 
an  accurate  RANS  computation,  but  that  transition  to  LES  models  away  from  the 
wall.  This  technique  alleviates  the  difficulties  of  prescribing  matching  conditions 
present  in  the  previous  methods.  The  most  well  known  approach  of  this  type  is  that 
of  Spalart  et  al.  (1997)  called  detached-eddy  simulation  (DES).  This  method  uses 
a  modified  formulation  of  the  one-equation  Spalart-Allmaras  (SA)  eddy  viscosity 
model.  Unlike  the  RANS  version,  this  model  uses  a  length  scale  that  is  the  distance 
from  the  wall  in  the  near-wall  region  and  switches  to  the  LES  filter  width  away  from 
the  wall.  DES  was  originally  conceived  for  massively  separated  flows  with  LES 
resolving  the  separated  region  while  RANS  computes  the  boundary  layer.  Further 
investigations  of  this  approach  examined  the  use  of  this  technique  in  plane  channel 
flow  of  varying  Reynolds  numbers  to  determine  how  it  would  behave  in  flows  without 
separation  (Nikitin  et  al.,  2000).  Some  encouraging  results  were  obtained,  as  the 
viscous  sub-layer  and  near- wall  logarithmic  profile  were  well  predicted  over  a  range 
of  Reynolds  numbers.  However,  the  skin  friction  coefficient  was  under-predicted 
by  about  15%  due  to  the  development  of  a  spurious  buffer  layer  in  the  logarithmic 
layer  which  shifted  the  mean  velocity  upwards.  Piomelli  et  al.  (2003)  were  able 
to  mitigate  this  problem  by  using  stochastic  forcing  (see  section  1.2.1  for  a  more 
complete  discussion).  A  similar  under-prediction  of  the  skin  friction  coefficient  was 
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observed  for  attached  flow  over  a  flat  plate  (Caruelle  and  Ducros,  2003).  In  addition, 
when  used  to  study  pressure  induced  separation  on  a  flat  plate,  DES  over-predicted 
the  length  the  separated  region  by  over  a  factor  of  2. 

While  the  DES  approach  has  been  extended  by  Strelets  (2001)  to  other  RANS 
turbulence  models,  it  is  also  possible  to  blend  viscosity  models  to  smoothly  transi¬ 
tion  from  RANS  near  the  wall  to  LES  away  from  it.  A  simple  example  from  Baggett 
(1998)  is 

Tij  ~  ^ TijSij  =  -  ((1  -  P(y))vLE S  +  P{y)v RANS)  Sij.  (1.1) 

In  this  equation,  the  blending  function  (3(y)  is  1  when  the  viscosity  is  purely  RANS 
and  0  when  it  is  purely  LES.  Typically,  f3  is  taken  to  be  1  at  the  wall,  followed  by 
a  smooth  transition  to  0  at  a  location  away  from  the  wall.  Above  this  location, 
the  simulation  uses  only  the  LES  viscosity.  While  this  approach  can  be  tuned  to 
yield  good  results  in  certain  situations,  (3  cannot  be  determined  theoretically  and  is 
expected  to  be  different  depending  on  the  numerical  method,  grid  resolution,  and 
SGS  models  used  in  a  given  computation.  We  are  also  unaware  of  any  technique  for 
dynamically  adjusting  /?.  Other  authors  (Germano,  1999;  Speziale,  1998;  Aruna- 
jatesan  and  Sinha,  2001)  have  also  worked  on  constructing  universal  models  that 
asymptotically  approach  RANS  or  LES  models  depending  on  the  grid  spacing  and 
flow  conditions,  all  with  limited  success. 

Some  authors  in  the  meteorological  community  use  a  different  technique  that  is 
similar  to  a  blended  eddy  viscosity  model.  An  extra  stress  is  added  to  the  Navier- 
Stokes  equations  with  a  prescribed  form  that  is  chosen  to  decrease  to  zero  at  some 
point  away  from  the  wall  (Brown  et  al,  2001).  This  gives  an  equation  for  the  stress 
to  be: 

ri2  =  -  J  Cca(y)  |u|  indy,  (1.2) 

with  a(y )  being  the  aforementioned  shape  function,  and  the  subscript  2  denoting 
the  wall-normal  direction.  The  equation  is  used  to  solve  for  the  index  i  =  1  (the 
streamwise  wall  stress)  and  i  =  3  (the  spanwise  wall  stress).  The  magnitude  of 
this  model  can  be  adjusted  with  Cc.  Cederwell  (2001)  chose  this  constant  to  match 
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experimental  observations  of  stress  in  a  tree  canopy,  and  it  can  also  be  tuned  to 
enforce  a  logarithmic  law  in  the  mean  velocity  profiles,  as  was  done  by  Chow  et  al. 
(2005).  The  difficulty  with  this  type  of  model  is  that  the  shape  must  be  adjusted 
by  trial  and  error,  and  in  the  meteorological  community,  these  functions  have  been 
adjusted  to  match  the  stresses  from  the  de  facto  rough  wall  present  in  environmental 
flows  consisting  of  trees,  rocks,  houses,  etc. 


1.2  Standard  Wall  Models 

After  examining  some  of  the  alternatives,  wall  models  will  now  be  considered.  Wall 
models  are  categorized  into  three  main  groups:  algebraic  models  that  use  a  simple 
relationship  between  the  wall  stress  and  LES  state,  two-layer  models  that  utilize 
some  set  of  near- wall  dynamics  to  prescribe  the  wall  stresses,  and  control-based 
wall  models  that  formulate  controllers  to  regulate  the  LES  via  wall  stress  inputs. 
In  the  notation  here,  a  standard  wall  model  will  denote  either  of  the  two  former 
approaches,  since  these  have  been  in  use  the  longest. 

An  additional  feature  these  models  share  is  that  they  aim  to  provide  boundary 
conditions  only  by  accounting  for  unresolved  physics.  This  is  typically  accomplished 
by  prescribing  wall  stresses  on  the  wall-parallel  velocity  components  while  the  wall- 
normal  velocity  is  set  to  zero.  This  restriction  arises  from  the  fact  that  it  is  difficult 
to  determine  an  appropriate  penetration  velocity  from  purely  physical  reasoning 
since  both  this  component  and  its  wall-normal  derivative  are  zero  at  the  wall.  An 
additional  difficulty  is  that  if  the  penetration  velocity  is  non-zero,  it  must  be  set 
such  that  there  is  no  net  mass  flux  through  the  wall.  This  means  that  it  will  not 
be  possible  to  determine  this  velocity  from  local  LES  data,  requiring  additional 
complexity  from  the  wall  model.  Therefore,  in  the  following  discussion  it  should 
be  understood  that  the  models  that  are  described  are  wall-stress  models  with  zero 
penetration  velocity.  This  will  not  be  the  case  when  control-based  wall  models  are 
discussed,  as  a  controller  can  provide  an  appropriate  penetration  velocity. 
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1.2.1  Algebraic  Wall  Models 

Wall-stress  models  were  the  first  type  of  near  wall  treatment  considered  for  LES. 
This  type  of  model  replaces  the  classical  no-slip  boundary  conditions  in  the  stream- 
wise  and  spanwise  directions  with  wall  stresses  so  that  the  near-wall  turbulence 
need  not  be  resolved.  The  first  attempt  at  such  a  model  was  by  Deardorff  (1970) 
who  used  the  following  model  in  a  LES  of  plane  channel  flow  at  infinite  Reynolds 
number: 


d2u  _  1  d2u 

dy 2  nyj  dz2 

d2w  d2w 
dy 2  dx 2 ’ 


(1.3) 

(1.4) 


with  u  and  w  being  the  filtered  streamwise  and  spanwise  velocity  components, 
respectively,  while  yi  is  the  location  of  the  first  grid  point  off  the  wall  and  k  is  the 
von  Karman  constant.  These  boundary  conditions  are  unique  in  that  they  impose 
a  condition  on  the  second  derivative  at  the  wall.  Note  that  in  the  mean,  these 
conditions  imply  a  logarithmic  profile  at  the  boundary.  When  combined  with  a 
no-penetration  condition  at  the  wall,  the  conditions  on  u  and  w  provide  all  the  wall 
data  required  by  the  simulation.  Using  this  model,  Deardorff  was  able  to  compute 
the  flow  in  a  plane  channel,  although  the  mean  statistics  were  not  in  good  agreement 
with  the  experimental  data.  This  deficiency  cannot  be  solely  attributed  to  the  wall 
model,  however,  as  the  grid  resolution  was  too  coarse  to  properly  resolve  even  the 
outer  length  scales. 

Schumann  (1975)  was  the  first  to  implement  what  is  now  considered  a  standard 
wall-stress  model  in  a  LES  of  plane  channel  flow.  The  wall  stresses  were  determined 
by  assuming  that  they  were  in  phase  with  the  velocities  at  the  first  interior  grid 
point,  and  that  the  local  deviation  from  the  mean  was  proportional  to  the  deviation 
from  the  mean  of  the  LES  velocity  at  the  nearest  wall-normal  grid  point. 
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Specifically,  the  following  model  was  used: 


w  ,  ,  \du 

TU  =  (»  +  V,)^ 


(r-) 
{«)  ivi) 


'32 


.  .  dw 


w 


(15) 

(1.6) 


with  (•)  denoting  plane  averaging,  v  is  the  molecular  viscosity,  and  vt  the  eddy  vis¬ 
cosity.  Also,  { tw )  represents  the  averaged  streamwise  wall  stress.  This  can  either 
be  taken  to  balance  the  applied  mean  pressure  gradient  (which  is  only  applicable 
in  channel  flow)  or  iteratively  solved  to  impose  that  the  plane-averaged  streamwise 
velocity  at  t/j,  the  first  grid  point  in  the  channel  interior,  satisfies  the  logarithmic 
law  of  the  wall  by  assuming  the  boundary  layer  is  in  equilibrium.  This  model 
produced  much  better  results  in  channel  flow  than  Deardorff’s  coarse  grained  cal¬ 
culations.  Several  improvements  have  been  suggested  to  this  type  of  model,  such  as 
the  method  by  Piomelli  et  al.  (1989)  (see  Section  1.2.3)  which  moved  the  matching 
point  downstream  to  account  for  the  inclination  of  near- wall  vortical  structures. 
Grotzbach  (1987)  used  a  model  of  this  type  to  impose  heat  fluxes  at  the  wall  in 
computations  involving  heat  transfer. 

As  mentioned  in  the  previous  section,  wall  modeling  has  also  been  of  great 
importance  in  simulating  environmental  flows  where  the  wall  stresses  are  typically 
set  based  on  enforcing  the  logarithmic  profile  locally  and  instantaneously  (Mason 
and  Callen,  1986).  Mean  velocity  profiles  other  than  the  logarithmic  law  have  also 
been  used  to  compute  the  wall  stress  in  (1.5).  The  work  of  Werner  and  Wengle 
(1991),  for  example,  used  a  near-wall  linear  profile  with  a  power  law  further  from 
the  wall.  The  predictions  of  these  computations  tend  to  be  similar  to  those  obtained 
by  Schumann  (1975)  and  Piomelli  et  al.  (1989). 

Mason  and  Thomson  (1992)  used  a  stochastic  backscatter  model  in  conjunction 
with  the  wall  model  of  Mason  and  Callen  (1986).  This  model  attempted  to  account 
for  the  effects  of  the  backscatter  of  energy  from  the  small  scales  to  the  large  scales 
by  adding  a  random  force  to  the  Navier-Stokes  equations  in  the  near-wall  region.  By 
adjusting  the  amplitude  of  this  force,  they  were  able  to  significantly  improve  upon 
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the  mean  velocity  profile  of  Mason  and  Callen  (1986).  Both  Mason  and  Thomson 
(1992)  and  Piomelli  et  al.  (2003)  reported  that  the  stochastic  force  “breaks  up”  the 
large  structures  and  produces  a  less  correlated  velocity  field.  The  exact  manner  in 
which  this  improves  the  prediction  of  the  mean  velocity  in  the  outer  layer  is  unclear, 
although  it  is  likely  that  the  random  forcing  adds  energy  which  is  transported  to 
the  outer  flow.  However,  it  is  clear  from  the  instantaneous  flow  contours  that  the 
resulting  flow  structures  do  not  correspond  to  the  well  known  features  in  attached 
boundary  layers.  In  addition,  there  is  currently  no  way  of  selecting  the  amplitude 
of  the  random  force  a  priori.  This  result  cannot  therefore  be  used  as  a  general 
purpose  wall  model,  but  does  provide  evidence  that  standard  wall  models  must  be 
corrected  in  order  for  a  good  prediction  of  the  mean  velocity  profile  to  be  obtained. 


1.2.2  Two-Layer  Wall  Models 


The  other  type  of  standard  wall  model  uses  simplified  versions  of  the  thin  boundary 
layer  equations  (TBLE)  to  determine  the  wall  stress.  These  equations  are  given  by: 


dui 

aF  + 


diiiUj 
d  x j 


1  dp  d  .  .  dui 

p  dxi  dy  dy 


(1.7) 


where  all  diffusion  terms  not  in  the  wall-normal  direction  are  assumed  to  be  small. 
The  boundary  conditions  for  (1.7)  are  taken  to  be  no-slip  at  the  wall  and  the  LES 
velocity,  ui7  at  the  outer  boundary  y  =  ym.  The  wall-normal  velocity  is  computed 
to  satisfy  the  continuity  equation 

,  .  fy  ( du  dw\  ,  .  . 

=  (fe+  &)<*/■  (1'8) 

for  y  e  (0,  ym]. 

By  neglecting  the  convective  terms  in  (1.7),  Hoffman  and  Benocci  (1995)  con¬ 
structed  a  local  model  by  integrating  the  TBLE  in  the  wall-normal  direction: 
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The  LES  pressure  was  used  in  these  equations  (assuming  that  pressure  is  constant 
in  the  wall- normal  direction).  The  time  derivative  was  evaluated  directly  from 
the  LES  computation  so  that  the  model  could  be  evaluated  instantaneously  and 
locally  without  having  to  store  the  TBLE  state.  Finally,  a  mixing  length  eddy 
viscosity  model  was  used  to  compensate  for  the  neglected  terms.  This  approach 
was  implemented  in  plane  and  rotating  channel  flow  with  reasonable  results. 

In  an  effort  to  incorporate  more  physics  into  the  wall-stress  models,  Balaras  et  al. 
(1994)  introduced  a  two-layer  approach  that  solves  an  additional  set  of  dynamical 
equations  near  the  wall.  The  near-wall  equations  are  solved  on  a  fine  wall-normal 
grid,  as  shown  in  Fig.  1.1.  The  wall  stress  computed  by  the  inner  layer  is  then  used 
as  a  boundary  condition  for  the  LES.  Balaras  et  al.  (1996)  attempted  a  model  of 
this  type  in  a  plane  channel,  square  duct,  and  rotating  channel  using  the  full  TBLE 
to  compute  a  near- wall  velocity  field  on  a  fine  mesh  embedded  in  the  first  cell  of 
the  LES  grid.  Savings  over  the  full  LES  equations  are  realized  since  the  TBLE  grid 
need  only  be  refined  in  the  wall-normal  direction  and  uses  the  LES  grid  spacing  in 
the  wall  parallel  directions.  Further,  since  the  pressure  is  applied  from  the  LES  and 
v  is  solved  to  satisfy  continuity,  no  pressure  solution  is  required  for  the  near-wall 
region.  Note  also  that  an  eddy  viscosity  model  is  often  used  to  compensate  for  the 
neglected  terms  and  the  large  wall-parallel  grid  spacing,  and  most  practitioners  use 
some  form  of  a  mixing  length  model  with  near  wall  damping. 

Since  this  was  the  same  approach  taken  by  Hoffman  and  Benocci  (1995),  but 
without  neglecting  non-linear  terms,  it  was  unclear  what  physics  the  near- wall 
model  should  retain  and  what  could  be  neglected.  Cabot  (1996)  considered  a  variety 
of  different  near-wall  models  in  an  LES  of  a  backward  facing  step  flow.  The  results 
were  mixed  for  each  model.  In  particular,  some  quantities,  such  as  the  pressure 
coefficient  after  the  step,  were  poorly  predicted  by  all  the  models.  More  recently, 
Wang  and  Moin  (2002)  used  a  two-layer  model  to  compute  an  airfoil  trailing  edge 
flow.  Several  variants  of  (1.7)  with  a  mixing-length  eddy  viscosity  were  considered: 
1)  setting  wall-normal  diffusion  equal  to  zero,  2)  wall- normal  diffusion  balancing  the 
LES  pressure  gradient,  and  3)  the  full  TBLE  equation.  In  case  3),  they  dynamically 
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LES  grid  (outer  scales) 


Approximate  B.C’s 
(l  w/>  tw3) 


Figure  1.1:  Two-layer  model  schematic  (courtesy  of  M.  Wang). 

adjusted  the  coefficient  of  the  mixing-length  eddy  viscosity  model  to  match  the  LES 
and  RANS  shear  stresses  at  the  interface.  The  most  accurate  results  were  achieved 
by  the  last  approach.  Mean  velocity  profiles  and  skin  friction  are  presented  in  Fig. 
1.2,  which  show  good  agreement  with  a  resolved  LES.  However,  when  Catalano 
et  al.  (2003)  used  case  2  to  compute  flow  over  a  cylinder  at  high  (super-critical) 
Reynolds  numbers,  the  Reynolds  number  dependence  of  the  drag  coefficient  was 
not  captured.  Problems  were  also  encountered  when  using  very  coarse  grids  in 
the  trailing  edge  simulation.  These  results  illustrate  the  primary  difficulty  with 
standard  wall  models.  Although  some  success  has  been  obtained  using  them  in 
simple  geometries  at  low  to  moderate  Reynolds  numbers,  none  has  demonstrated 
the  robustness  needed  to  be  used  on  very  coarse  grids  at  high  Reynolds  numbers. 

The  one-dimensional  turbulence  (ODT)  model  of  Kerstein  et  al.  (2001)  was 
recently  used  as  a  SGS  and  wall  model  for  pressure-driven  plane  channel  flow  by 
Schmidt  et  al.  (2003).  In  order  to  apply  this  model  to  fully  three-dimensional  flow, 
the  standard  model,  which  only  includes  wall-normal  diffusion,  was  augmented  to 
include  the  LES  pressure  gradient  and  a  convection  term  similar  to  (1.7),  only  here 
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Figure  1.2:  Trailing-edge  skin  friction  coefficient  (left)  and  mean  velocity  profiles 

(right); - :  dynamic  k, - :  constant  k  —  0.4  (left)  and  resolved  LES  (right), 

.  :  resolved  LES  (left),  •:  experiment  of  Blake  (1975). 


the  convecting  velocity  is  taken  from  an  average  over  the  LES  time  step  in  the 
cells  in  which  the  ODT  model  is  used.  The  ODT  is  advanced  using  a  smaller  time 
step  to  include  “eddy  events”:  random  perturbations  to  the  velocity  designed  to 
mimic  turbulent  eddies.  To  couple  the  ODT  to  the  LES,  these  events  were  allowed 
to  extend  out  into  the  LES  over  its  first  few  grid  points.  Reasonable  results  were 
reported  over  a  variety  of  Reynolds  numbers,  although  the  slope  of  the  logarithmic 
profile  becomes  increasingly  over-predicted  with  increasing  Reynolds  number.  Also, 
some  discrepancies  are  noted  in  the  wake  region  of  the  flow,  particularly  in  the  wall- 
normal  rms  velocity  fluctuations.  An  additional  issue  is  the  high  computational 
expense  of  the  method. 

1.2.3  Deficiencies  of  Standard  Wall  Models  in  High  Reynolds 
Number  Flow 

The  previous  section  illustrated  that  many  variants  of  wall  stress  models  have  been 
proposed  over  the  past  thirty-five  years.  In  plane  channel  flow,  all  of  these  models 
provide  streamwise  and  spanwise  wall  stresses  at  each  grid  point  on  the  wall  while 
retaining  a  zero  penetration  velocity.  In  most  cases,  either  an  averaged  or  instan¬ 
taneous  logarithmic  profile  is  used  to  predict  the  mean  wall  stress.  Before  moving 
on  to  control-based  wall  modeling  techniques,  it  is  useful  to  evaluate  these  models 
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in  the  test  case  that  will  be  considered  in  this  work. 

The  first  model  considered  will  be  the  shifted  model  of  Piomelli  et  al.  (1989). 
If  we  denote  the  streamwise  wall  stress  by  r™2  and  the  spanwise  wall  stress  by  r^, 
the  shifted  model  specifies  the  local  stresses  by 


r^Or,  z)  = 


u(x  +  S,  yuz ) 
(u) 

w(x  +  6,  yi,  z) 

(«> 


(rw) 
< rw ) 


(1.10) 

(1.11) 


with  (rw)  being  the  mean  streamwise  wall  stress  computed  by  assuming  a  log- 
law  velocity  profile  near  the  wall.  Recall  that  in  this  and  all  other  models,  the 
transpiration  velocity  is  taken  to  be  zero. 

The  other  model  examined  is  a  variant  of  the  TBLE  equation  model  (1.7).  The 
following  equations,  as  presented  by  Wang  and  Moin  (2002),  are  used: 

^  +  ^=0,  <  =  1.3  (1.12) 

with  the  mixing-length  eddy  viscosity  model 

^  =  ny+  (l  -  e~y+/Ay  (1.13) 


with  k  =  .41  and  A  =  17.9.  In  this  form,  the  TBLE  model  is  simply  an  ODE  that 
can  be  analytically  integrated  to  yield 

<v)  =  c'[t^, W)iy'+C2-  (U4) 

The  integration  constants  are  set  such  that  C2  =  0,  enforcing  the  no-slip  condition, 
and  Ci  is  found  by  the  matching  condition  Uj(yi)  =  w^les^i)-  The  wall  stress  is 
directly  identifiable  as  C\. 

Figure  1.3  demonstrates  that  both  models  perform  nearly  identically  in  channel 
flow  at  ReT  =  4000  using  a  uniform  mesh  of  32  x  33  x  32  grid  points  and  the  agree¬ 
ment  with  the  standard  logarithmic  law  is  not  satisfactory.  Additional  evidence  of 
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this  insensitivity  is  offered  by  considering  an  extremely  simple  wall  model: 

Ta  =  puT{ui(yx)  -Ui)  +  h^t  (1.15) 

for  i  =  1,3,  where  the  friction  velocity  is  defined  by  <  >=  pui  and  h  is  the 

channel  half  height.  In  this  model,  Ui  is  a  matching  velocity  set  a  priori ,  in  this  case 
to  match  the  logarithmic  profile  at  yi,  and  the  second  term  on  the  right-hand  side 
is  present  to  balance  the  mean  pressure  gradient.  This  model  can  be  seen  to  be  a 
simple  feedback  control  setting  the  wall  stress  to  target  a  mean  value  for  the  velocity. 
In  fact,  it  is  even  simpler  than  a  typical  feedback  controller  since  the  gain  is  naively 
taken  to  be  unity.  However,  the  mean  profile  it  produces  when  U\  matches  the  law 
of  the  wall  and  =  0,  as  shown  in  Fig,  1.3,  is  almost  identical  to  the  other  two 
models  that  use  advanced  techniques  and  knowledge  of  turbulent  flows  to  predict 
the  wall  stress.  It  is  reasonable  to  suppose  that,  despite  their  differences,  all  the 
models  have  an  underlying  structure  that  give  the  same  wall  stress  predictions.  It 
seems  clear  that  a  model  based  on  these  principles  will  encounter  difficulties  in  flows 
at  high  Reynolds  numbers  on  coarse  grids. 

1.3  Control-Based  Wall  Models 

Algebraic  and  TBLE  wall  models  have  produced  successes  in  certain  cases,  but 
none  has  been  demonstrated  to  be  robust  enough  to  be  used  in  a  general  setting. 
This  is  likely  due  to  standard  models  relying  on  obtaining  the  missing  physics 
from  coarse  simulations  without  addressing  the  effects  of  SGS  modeling  errors  and 
numerical  errors  present  near  the  wall.  Cabot  (1997)  provided  direct  evidence  of  the 
significance  of  these  errors  by  using  the  wall-stresses  obtained  from  a  resolved  LES 
of  a  backward  facing  step  as  a  wall  model.  These  were  then  used  in  an  LES  with 
the  same  initial  conditions  and  resolution  away  from  the  wall,  but  with  the  first 
ten  near- wall  points  removed.  The  results  demonstrated  that  even  the  “correct” 
wall  stresses  could  not  produce  a  wall  model  that  was  more  accurate  than  standard 
phenomenologically  derived  techniques.  What  was  needed  was  a  method  that  could 
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Figure  1.3:  Mean  flow  profile  using  standard  wall  models; .  :  shifted  model  of 

Piomelli  et  al.  (1989), - :  algebraic  model  of  Wang  and  Moin  (2002), - : 

simple  wall  model  (1.15), - :  logarithmic  profile  ( u+  =  2.41  log y+  +  5.2). 
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actively  regulate  an  LES.  The  first  attempts  at  such  a  technique  were  by  Nicoud 
et  al  (2001),  who  in  fact  tried  two  different  approaches. 

The  first  approach  involved  the  application  of  optimal  control  theory  to  imple¬ 
ment  a  regulator  to  provide  the  wall  stresses,  since  it  is  unknown  how  to  compensate 
for  numerical  and  SGS  errors.  The  approach  used  was  similar  to  that  of  optimal 
flow  control  (Bewley  and  Moin,  1997).  A  cost  function  was  defined  that  measured 
the  plane-averaged  deviation  of  the  LES  velocity  from  that  of  the  logarithmic  pro¬ 
file.  Adjoint  equations  (see  Chapter  2)  were  used  to  compute  the  gradient  of  this 
cost  function  with  respect  to  the  control,  in  this  case  taken  to  be  the  streamwise 
and  spanwise  wall  stresses  (the  transpiration  velocity  was  set  to  zero).  Several 
approximations  were  made  in  the  formulation  of  the  adjoint  equations,  as  well 
as  the  LES-based  equations  used  to  compute  the  physical  state  required  for  the 
adjoint’s  solution.  Further,  the  controls  were  optimized  over  one  time  step  only. 
These  reductions  imply  that,  while  optimal  control  theory  was  used,  the  controller 
was  in  fact  sub-optimal.  Despite  this,  the  results  from  a  LES  of  channel  flow  at 
ReT  =  4000  using  32  x  33  x  32  cells  display  a  good  prediction  of  the  mean  velocity 
profile  throughout  the  domain  (Fig.  1.4). 

The  second  regulator  implemented  by  Nicoud  et  al.  was  a  feedback  controller 
constructed  from  the  results  of  the  sub-optimal  control.  A  linear  stochastic  esti¬ 
mation  (LSE)  (Bagwell  et  al,  1993)  was  performed  on  the  wall  stresses  produced 
by  the  sub-optimal  control  framework  to  determine  the  optimal  linear  correlation 
between  the  velocity  field  and  the  wall  stresses.  The  resulting  controller  was  then 
of  the  form  of  a  kernel  convolved  with  the  velocity  field.  Results  of  this  regulator 
at  ReT  =  640  and  20  000  showed  a  good  prediction  of  the  mean  velocity  profile. 

Both  the  LSE  and  the  suboptimal  regulators  were  extended  by  Baggett  et  al 
(2000).  First,  transpiration  velocity  was  added  to  the  control  set  of  the  sub-optimal 
regulator.  However,  this  addition  did  not  significantly  improve  the  model’s  predic¬ 
tions  relative  to  the  improvement  obtained  when  replacing  a  standard  wall  model 
with  a  control-based  wall  model.  A  cost  function  including  terms  measuring  the 
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Figure  1.4:  Mean  flow  profile  using  a  control-based  wall  model; . :  shifted  model 

of  Piomelli  et  al.  (1989), - :  control-based  model  of  Nicoud  et  al.  (2001), - : 

logarithmic  profile. 
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deviation  of  the  rms  velocity  fluctuations  were  also  considered.  The  rms  target  pro¬ 
files  were  taken  from  the  LES  of  Kravchenko  et  al.  (1996)  using  zonally  embedded 
meshes.  While  minimization  of  this  cost  function  did  slightly  improve  the  match 
between  the  predicted  rms  velocity  fluctuations  and  those  of  Kravchenko  et  al. ,  a 
decrease  in  the  accuracy  of  the  prediction  of  the  mean  velocity  also  occured.  This 
is  possibly  due  to  the  control  objectives  being  in  conflict  with  each  other.  Baggett 
et  al.  also  further  investigated  the  use  of  the  LSE  feedback  regulator  by  using  the 
one. previously  obtained  by  Nicoud  et  al.  in  new  channel  flow  simulations  with 
different  numerical  methods.  When  implemented  in  a  code  using  fourth  order  fi¬ 
nite  differencing  to  evaluate  the  spatial  derivatives  (Nicoud  et  al.  used  a  second 
order  formulation),  the  mean  profile  was  not  as  well  predicted.  This  indicates  the 
controller  was  adjusting  the  wall  stresses  based  on  the  discretization  stencil  used 
in  the  simulation.  An  even  greater  change  was  observed  when  the  SGS  model  of 
Cabot  and  Moin  (2000)  was  used  to  increase  the  eddy  viscosity  in  the  first  cell.  In 
this  case,  the  slope  near  the  wall  was  significantly  over-predicted,  resulting  in  the 
intercept  of  the  law  of  the  wall  being  too  great.  This  result  shows  that  the  con¬ 
troller  strongly  reacts  to  the  SGS  model.  Specifically,  it  increases  the  rms  velocities 
to  compensate  for  the  SGS  model  not  carrying  enough  turbulent  stress. 

The  results  of  Nicoud  et  al.  (2001)  and  Baggett  et  al.  (2000)  indicate  that  each 
simulation  will  require  its  own  active  controller,  at  least  for  a  sufficient  time  to 
derive  an  LSE-based  regulator  which  can  compensate  for  the  numerical  and  SGS 
modeling  errors  present  in  the  simulation.  In  addition  to  the  significant  cost  of  the 
sub-optimal  control,  its  extension  to  more  complex  flows  is  limited  by  the  need  to 
have  a  target  mean  velocity  profile  known  a  priori ,  i.e.  the  method  is  not  predictive. 

These  issues  led  Templeton  et  al.  (2002)  to  propose  a  different  type  of  feedback 
regulator.  This  approach  uses  a  near-wall  model  similar  to  Wang  and  Moin  to 
generate  target  velocity  profiles.  Since  these  models  are  valid  only  near  the  wall,  the 
cost  function  is  similarly  only  defined  in  this  region.  To  reduce  the  computational 
cost,  a  predictor-corrector  approach  was  used  in  that  the  wall  stress  determined 
by  the  near-wall  model  was  used  as  an  initial  guess  for  the  control.  Then,  one 
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optimization  iteration  was  performed  by  descending  along  the  gradient  direction 
of  the  cost  function.  Since  an  adjoint  equation  is  difficult  to  formulate  for  the 
trailing  edge  airfoil  flow,  significant  approximations  were  made  to  the  gradient 
such  that  it  was  computed  using  surface  data  only  (Mohammadi  et  al. ,  2000). 
Unfortunately,  this  approximation  proved  too  severe,  indicating  that  a  significant 
amount  of  accuracy  in  the  gradient  is  required  for  a  successful  regulator. 


1.4  Research  Objectives 

There  are  two  outstanding  issues  in  the  development  of  active  LES  regulators:  pre¬ 
dictability  and  cost.  The  latter  is  a  significant  issue  because  the  purpose  of  a  wall 
model  is  to  reduce  computational  expense  to  make  simulating  high  Reynolds  num¬ 
ber  flows  more  tractable.  In  the  work  involving  sub-optimal  control,  the  cost  of 
the  wall  model  is  on  the  order  of  ten  times  the  cost  of  the  rest  of  the  simulation. 
This  occurs  because  both  the  adjoint  and  LES  equations  must  be  solved  once  per 
iteration,  and  0(10)  iterations  are  required  to  obtain  a  converged  solution.  There¬ 
fore,  one  objective  of  this  work  was  to  reduce  the  computational  effort  required  per 
iteration  to  enable  the  model  to  be  used  efficiently. 

The  other  issue  that  must  be  resolved  is  the  predictability  of  the  method.  In  the 
work  of  Nicoud  et  al. ,  the  target  profile  used  was  prescribed  a  priori  While  this 
can  be  done  in  canonical  boundary  layers  since  the  mean  velocity  is  known,  in  an 
arbitrary  flow  the  mean  velocity  profile  will  not  be  known  before  the  computation 
is  performed.  The  problem  of  predictability  will  be  addressed  through  the  use  of 
RANS  equations  to  determine  the  target  profile. 

Chapter  2  will  derive  the  continuous  formulation  of  all  the  equations  needed 
in  this  work.  Issues  related  to  the  numerical  solution  of  these  equations  will  be 
presented  in  Chapter  3,  with  special  emphasis  on  techniques  to  discretize  the  adjoint 
equations  and  the  choice  of  cost  functions  consistent  with  those  discretizations. 
This  will  be  followed  by  Chapter  4  in  which  an  efficient  method  for  solving  the 
optimization  problem  in  plane  channel  flow  will  be  presented.  In  order  to  make 
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this  approach  predictive,  Chapter  5  will  demonstrate  how  RANS  velocity  profiles 
can  be  incorporated  into  the  cost  function  definition.  Some  final  thoughts  and 
conclusions  will  be  offered  in  Chapter  6. 

1.5  Accomplishments 

•  Evaluated  the  applicability  of  cost  function  gradients  computed  using  the 
method  of  incomplete  sensitivities  to  the  problem  of  control-based  wall  mod¬ 
eling  (see  Templeton  et  al.  (2002)). 

•  Demonstrated  the  importance  of  the  interplay  between  cost  function  definition 
and  adjoint  discretization  in  constructing  an  accurate  sub-optimal  regulator 
(Chapter  3). 

•  Significantly  reduced  the  computational  expense  of  the  optimal  control-based 
wall  model  by  taking  advantage  of  the  adjoint  formulation  and  cost  function 
structure  (Chapter  4). 

•  Determined  a  method  to  incorporate  RANS  velocity  profiles  into  the  cost 
function  definition  to  make  wall  models  based  on  optimal  control  techniques 
predictive  (Chapter  5). 


Chapter  2 


Derivation  of  the  Continuous 
Equations 


2.1  Introduction  and  Notation 

In  this  chapter  we  present  the  continuous  equations  that  will  be  considered  in 
this  work.  The  first  set  of  equations  presented  will  be  the  incompressible  Navier- 
Stokes  equations.  For  convenience,  we  will  define  q  =  [u,  p]  to  represent  the  full 
state.  In  what  follows,  the  velocity  u  will  be  written  interchangeably  as  (u\,U2,  U3) 
or  ( u,v,w ),  which  represent  the  components  of  the  velocity  field  in  the  (x\,X2,  £3) 
directions,  respectively.  It  will  often  be  convenient  to  refer  to  the  coordinate  axes  as 
(x,y,z).  When  considering  velocity  components  individually,  the  notation  ( x,y,z ) 
and  (u,v,w)  will  be  utilized.  In  this  work,  summation  over  repeated  indices  (i.e. 
i,j,  etc.)  is  implied,  except  when  specifically  indicated. 

The  first  set  of  equations  that  will  be  presented  are  those  for  incompressible, 
Newtonian  fluid  flow  with  constant  density.  Next,  the  LES  equations  are  derived, 
which  retain  the  large  scales  of  the  flow  while  modeling  the  small  ones.  Phys¬ 
ical  boundary  conditions  can  be  prescribed  for  the  Navier-Stokes  equations,  but 
transferring  these  conditions  to  the  LES  equations  can  present  some  computational 
difficulties. 
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The  goal  of  this  work  is  to  use  optimal  control  techniques  to  remedy  these  dif¬ 
ficulties  by  formulating  a  wall  model,  or  alternate  set  of  LES  boundary  conditions. 
In  order  to  use  such  techniques  efficiently,  the  adjoint  equations  of  the  LES  system 
must  be  derived.  The  solution  of  these  equations  can  be  thought  of  as  representing 
the  sensitivities  of  the  flow  to  disturbances,  or  of  being  Lagrange  multipliers  that 
account  for  the  constraint  of  the  LES  system  on  the  optimization  process.  In  the 
process  of  constructing  these  equations,  the  LES  equations  will  be  formally  lin¬ 
earized.  The  solution  of  the  linearized  equations  will  be  denoted  by  q'  =  \u' ,  p'\, 
where  v!i  corresponds  to  the  linearized  state  associated  with  u*.  Similar  notation 
will  be  used  for  the  adjoint  state,  q*  =  [u* ,  p*],  where  each  physical  variable  will 
have  a  corresponding  adjoint  variable. 

Since  much  of  this  work  involves  the  use  of  adjoint  equations  and  optimal  con¬ 
trol  techniques,  it  will  be  beneficial  to  express  many  of  the  equations  in  operator 
notation.  In  all  cases,  a  non-linear  operator  acting  on  a  vector  will  be  written  as 
A{q),  while  a  linear  operator  will  be  denoted  as  Bq. 


2.2  Navier-Stokes  Equations 

The  Navier-Stokes  operator,  which  is  used  to  write  the  equations  that  govern  in¬ 
compressible,  Newtonian  flows,  can  be  written  as: 

dpuj  ,  dpUjUj  .  dp  _JL.ii  ( 2Hi  — i\ 

dt  dxj  dxi  dxj  •  y  dxj  dx{  J 

duj 
dxj 

This  operator  defines  the  differential  operations  that  are  applied  to  the  state  q.  The 
Navier-Stokes  equations  can  be  written  compactly  as 

Af(q)  =  f,  (2.2) 

— # 

where  /,  a  vector  with  four  entries  at  each  spatial  and  temporal  location,  is  the 
source  term.  The  fourth  entry,  corresponding  to  the  divergence  operator  in  (2.1), 


(2.1) 
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must  be  everywhere  zero  to  enforce  the  divergence-free  constraint  on  the  velocity 
field.  The  other  terms  in  /  represent  momentum  sources,  which  can  come  from  the 
physics  or  be  control  inputs  into  the  system. 

The  final  component  required  to  define  the  Navier-Stokes  system  are  initial  and 
boundary  conditions.  The  velocity  and  pressure  fields  are  considered  to  exist  on  the 
closed  set  Cl,  while  the  Navier-Stokes  equations  are  valid  on  the  open  set  Cl  C  R3. 
The  boundary  of  the  set  is  defined  as 

dCl  =  ti\Cl. 

Without  loss  of  generality,  the  system  can  be  taken  to  start  at  t  =  0,  and  hence 
the  temporal  domain  is  (0,  T].  Therefore,  the  initial  and  boundary  conditions  are 
defined  as 


g|t«o  -  ?o(f) 
g(t,  x,  q :  x  €  3£2)  =  0. 


(2.3) 

(2.4) 


Note  that,  similarly  to  the  source  term  /,  g  can  also  contain  control  inputs  to  the 
system. 

We  denote  the  dimensionless  value  of  quantity  a  by  a*,  and  so  each  dimensionless 
variable  is  defined  by: 


u 

—  =  u' 
U 


-  ,,t 


P 

pU2 


x 


-  ~t 


—  =  x’ 
D 


D  ’ 


where  U  is  the  chosen  velocity  scale  and  D  the  chosen  length  scale.  Substituting 
these  expressions  into  (2.1)  yields  the  dimensionless  Navier-Stokes  operator: 


{du\  duj 
dx 

5u| 


d  X 
dxl  Re 


M<W)  = 


(2.5) 
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where  there  is  now  only  one  dimensionless  parameter,  the  Reynolds  number, 

Re=^. 

To  ease  the  notation,  all  quantities  should  be  taken  to  be  dimensionless  (without 
special  designation)  unless  otherwise  noted.  The  scales  used  to  make  the  variables 
dimensionless  will  be  presented  as  they  appear. 


2.3  Large-eddy  Simulation 


To  construct  the  LES  equations,  a  low-pass  filter  is  applied  to  the  state  q  in  order 
to  remove  the  small  scales.  The  filtered  q  is  denoted  by  q.  The  effects  of  the  small 
scales  on  the  large  ones  must  be  modeled.  The  LES  operator  is  written  as: 


m  = 


duiUj 

dxj 


■  dp _ d_  J_ 

dxi  dxj  Re 


(§ui  ,  duA  ,  drij_ 

y  dxj  '  Ox,  J  '  dxj 


where  is  called  the  sub-grid  scale  (SGS)  stress  and  is  given  by 


(2.6) 


Tij  =  UiUj  —  UiUj ,  (2.7) 

which  must  be  modeled  based  on  the  LES  state. 

In  this  work  we  will  use  an  eddy  viscosity  models  for  r^: 

=  2  vtSij  (2.8) 

with  §ij  being  the  filtered  strain  rate  tensor, 

q  _  1  (dui  , 
ij  2  \dxj  dxj  ’ 

and  vt  the  SGS  eddy  viscosity. 

A  common  SGS  eddy  viscosity  model  is  the  Smagorinsky  eddy  viscosity  model 
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(Smagorinsky,  1963): 

«  =  CSA 2 15| , 

where  Cs  is  a  model  coefficient,  A  is  the  filter  width,  and 


|s|  =  v4s,.s,. 


(2.9) 


The  Dynamic  model  (Germano  et  al,  1991;  Lilly,  1992)  allows  Cs  to  be  computed 
from  the  resolved  velocity  field: 


[MklMkly 


(2.10) 


where 


Mtj  =  A2|S|  Sij  -  A2f 

Lij  —  UiUj  lijUj 


(2.11) 

(2.12) 


and  [•]  is  an  averaging  operator.  In  (2.11),  :  denotes  a  test  filter  with  filter  width 
Ap  >  A.  In  flows  with  homogenous  directions,  the  averaging  operator  can  be 
applied  over  these  directions.  If  this  is  not  the  case,  the  dynamic  localization 
procedure  of  Ghosal  et  al.  (1995)  can  be  used  to  compute  the  model  coefficient. 
This  model  has  been  successfully  tested  in  a  range  of  applications  and  requires  no 
parameters  that  are  set  a  priori. 


2.4  Derivation  of  the  Adjoint  Operator 

In  this  section,  the  adjoint  operator  will  be  derived  from  the  LES  equations.  This  is 
in  contrast  to  Nicoud  et  al.  (2001)  in  which  the  adjoint  equations  were  formulated 
after  the  state  equation  was  already  discretized  in  time.  The  temporal  discretiza¬ 
tion  used  in  that  analysis  was  not  consistent  with  the  discretization  actually  used 
to  advance  the  state  equation.  In  contrast,  we  wish  to  determine  what  temporal 
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discretization  or  discretizations  will  perform  the  best  in  this  context,  with  perfor¬ 
mance  being  measured  by  the  accuracy  of  the  flow  statistics,  evaluation  effort,  and 
amenability  to  approximation  (to  further  reduce  the  evaluation  effort).  However, 
in  order  to  consider  the  temporal  discretization  of  the  adjoint  equations,  it  is  first 
necessary  to  derive  them  in  continuous  form.  This  has  been  done  many  times  be¬ 
fore  for  slightly  different  formulations  of  the  state  equations  (see  e.g.  Bewley  et  al. 
(2001)),  but  here  we  do  it  from  (2.6). 

Note  that  one  could  start  with  the  fully  discretized  Navier-Stokes  equations 
as  solved  in  the  code  and  then  derive  the  adjoint  equations  based  solely  on  them 
without  ambiguity.  While  this  would  achieve  the  best  possible  accuracy  because 
the  exact  gradient  of  the  discrete  system  could  be  computed,  the  evaluation  effort 
would  be  great  and  potential  for  approximation  would  be  quite  poor  (Nadarajah  and 
Jameson,  2000).  Further,  the  effort  required  to  construct  them  must  be  repeated 
for  each  new  code,  making  it  difficult  to  develop  a  general  wall-model  subroutine. 
In  this  section,  only  the  adjoint  operator  will  be  determined.  This  is  because  the 
adjoint  equations  themselves  will  be  determined  by  the  specific  control  problem 
under  investigation.  Thus,  one  should  look  at  the  adjoint  operator  as  a  system 
that  measures  the  sensitivity  of  the  flow  to  an  arbitrary  perturbation,  and  adjoint 
equations  as  governing  the  response  of  the  system  to  a  specific  perturbation. 

In  order  to  find  the  adjoint  equations,  it  is  necessary  to  first  write  the  equation 
for  the  Frechet  derivative  of  the  state  with  respect  to  a  small  perturbation  to  the 
control  input,  <f>.  The  Frechet  derivative  of  an  arbitrary  function,  J,  is  defined  to 
be 


DJ 7  M  + 
D(f>  e— *o  e 


where  (p  is  an  arbitrary  test  function.  Here,  the  gradient  of  J  with  respect  to  cp, 
DJ/ D<p,  is  a  linear  operator  acting  on  the  perturbation  (p.  See  Luenberger  (1969) 
for  additional  mathematical  details. 

The  work  here  will  focus  on  taking  the  Frechet  derivative  of  the  LES  state 
governed  by  operator  (2.6)  with  respect  to  a  control  input.  This  input  could  come 
through  the  initial  conditions,  boundary  conditions,  or  source  term.  To  ease  the 
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notation,  we  introduce  the  shorthand: 


(2.13) 


Taking  the  Frechet  derivative  of  Af(q)  yields  the  linearized  Navier-Stokes  operator: 


w=< 


du'- 

_ i 

dxj 


)  (£j  +  it)) 


(2.14) 


where  Al'q  is  a  linear  state  equation  acting  on  the  linearized  state  q'  about  a  base 
state  of  q ,  and  v  =  1/Re.  Here  we  have  ignored  the  sensitivity  of  vt  to  changes 
in  (f>.  This  approximation  was  shown  to  be  reasonable  for  short  time  intervals  by 
Chang  and  Collis  (1999). 


Frechet  differentiation  can  also  be  applied  to  (2.3)  and  (2.4)  to  identify  the 
initial  and  boundary  conditions  for  the  linearized  system: 


q'(t  =  0,x  eQ)  =  q'0  (2.15) 

g'(t,x,q'  :  x  G  dfi)  —  0.  (2-16) 

Note  in  (2.16)  g'  represents  the  boundary  conditions  for  the  linearized  system.  To 
write  the  linearized  Navier-Stokes  equations,  it  only  remains  to  take  the  Frechet 
derivative  of  the  LES  equations  to  obtain 

=  (2.17) 

As  was  mentioned  earlier,  the  specific  adjoint  equations  cannot  be  found  until  the 
optimization  problem  is  stated,  in  contrast  to  the  linearized  equations  which  are 
fully  known  once  the  LES  equations  are  prescribed. 

The  next  step  in  developing  the  adjoint  equations  is  to  determine  the  inner 
product  that  defines  the  space  in  which  the  functions,  q1,  exist.  Therefore,  we  take 
each  element  of  q'  to  be  a  function  in  £2(0  x  [0,T]).  The  inner  product  on  this 
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space  of  o,6  €  x  [0,T])  is  then 

(a,  6)  =  f  f  a(x,t)b(x,t)  dxdt.  (218) 

Jo  Jn 

Another  way  of  identifying  these  vectors  is  to  state  that  6  is  in  the  space  of  bounded 
linear  functionals  of  x  [0, T]),  denoted  by  /^(Q  x  [0,  T]),  which  is  the  dual  of 
the  original  space.  In  this  particular  case,  the  dual  and  original  spaces  are  the  same, 
and  so  6  is  also  an  element  of  the  original  space.  For  more  information  concerning 
dual  spaces  and  the  role  they  play  in  optimization,  the  interested  reader  is  referred 
to  Luenberger  (1969). 

In  this  formulation,  the  state  q'  is  a  member  of  the  original  space,  while  the 
adjoint  state,  q*,  is  as  yet  an  undetermined  element  in  the  dual  space.  The  adjoint 
operator  is  then  the  linear  operator,  Af* ,  acting  on  q* ,  that  satisfies  the  following 
identity 

=  +  BT,  (2.19) 

with  BT  being  terms  that  are  only  evaluated  at  the  temporal  and  spatial  boundaries 
of  the  domain.  Such  an  operator  is  guaranteed  to  exist  (Luenberger,  1969).  In  the 
case  of  differential  operators,  Gauss’  theorem  is  used  to  move  the  partial  derivatives 
from  q'  to  q*,  which  results  in  the  addition  of  the  boundary  terms,  BT,  which  will 
be  discussed  following  the  presentation  of  the  adjoint  operator. 
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Deriving  the  adjoint  equations  is  an  exercise  in  integration  by  parts  followed 
by  identifying  which  terms  are  in  which  adjoint  equation  by  comparing  them  with 
the  inner  product.  The  integration  by  parts  for  the  u'Af'qu  will  be  given  below  for 
clarity  (it  has  been  applied  twice  to  the  diffusion  terms): 


f  [  u*  ( +  u'u)  +  —  (uv'  +  u'v)  +  —  (mu'  +  u'w)  + 
Jo  J  q  \dt  dx  dy  dz 

d  ( du'  du'\  d  ( dv!  dv'\  d  ( dv!  dw'\\ 

—  -^~UT  (  "X b  -Q—  1  -  -X —  I  -  —  Vr  I  — - b  — —  I  > 

ox  \dx  ox  J  dy  \  oy 


dtf_ 

dx 


dx  J  dz  \dz  dx  J  J 


dxdt 


fT  f  (  ,du*  ,  ,  .du*  .  ,  , 

/  /  —  u  — - (uu  +  uu)- - - (vu  +  uv 

Jo  Jn\  dt  dx 


,d 

U—~UT 

dx 


du*  du* 

+ 


dx  dx 


.du*  .  .  ,  .du*  ,du* 

.  d  du*  .  d  du*  .  d  du* 

dy  dy  dx  dy  dz  dz 

d  du*  \ 

—  w'—vT-——  J  dxdt,  (2.20) 
dx  dz  J 


where  the  shorthand  uT  =  u  +  vt  has  been  used  to  write  the  equation  more  com¬ 
pactly.  The  notation  on  the  RHS  has  been  chosen  to  suggest  the  next  step  in 
deriving  the  adjoint  equations,  that  being  identifying  all  the  terms  multiplied  by  u! 
(including  those  coming  from  equations  multiplied  by  v* ,  w*  and  p*  terms  which 
were  not  shown)  into  the  equation  for  u*.  This  approach  is  implied  by  the  inner 
product  formulation.  Grouping  such  terms  yields  the  adjoint  operator: 

KV  =  I 


at 

duj 
dx i 


—  U 


du * 
'i  dxj 


—  U 


du’ 
j  dx. 


_  9e1 

dxj 


ET,  ((1/  +  1/<)  (sj+  s^)) 


.  (2.21) 


It  is  important  to  note  that  these  are  the  adjoint  equations  that  would  be  found 
following  Bewley  et  al.  (2001).  Nicoud  et  al.  (2001)  arrive  at  a  slightly  different  form 
by  applying  the  divergence-free  constraint  to  the  convective  terms  of  the  linearized 
system.  This  results  in  the  third  term  in  (2.21)  being 
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The  equivalence  between  the  two  formulations  is  found  by  taking  into  account  the 
adjoint  pressure  equation.  We  can  write: 

du*  du-jU*  duj 

The  resulting  source  term  in  the  adjoint  pressure  equation  is 

d  du*  _  d2UjU*  d  „  duj 
dxi d%i  dxf  dxi dx% 

Since  in  both  formulations  the  equation  for  p*  is  an  invertible  Laplace  equation,  we 
have 

P*i=Pl-  UjUj, 

where  p\  is  the  adjoint  pressure  following  Bewley  et  al.  and  p*2  is  the  same  quantity 
from  Nicoud  et  al..  When  this  is  used  in  (2.21),  it  is  seen  to  be  identical  to  the 
formulation  of  Nicoud  et  al.  (2001).  However,  their  formulation  requires  knowledge 
of  the  adjoint  source  term,  /*,  and  so  is  slightly  less  general.  Therefore,  the  formu¬ 
lation  presented  here  will  be  used,  although  if  the  divergence  of  the  adjoint  state  is 
constrained  to  be  zero,  the  two  approaches  are  equivalent. 

Of  course,  the  application  of  Gauss’  theorem,  which  resulted  in  the  derivatives 
being  moved  from  q'  to  q* ,  yields  a  series  of  terms  which  must  be  integrated  over 
the  spatial  and  temporal  boundaries.  To  specify  these  terms,  it  will  be  helpful 
to  introduce  some  new  notation.  First,  the  vector  nXi  is  the  component  of  the 
outward  facing  normal  vector  to  dQ,  in  the  xt  direction.  Second,  to  help  represent 
the  integration  by  parts  in  time,  denote 
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The  boundary  terms  arising  from  (2.20)  can  then  be  written  as 


BT„  =  u'u*  |o  +  f  (  u*  (( uvl '  +  uu')nx  +  (uv'  +  u'v)ny  +  (mu'  +  u'w)nz  +  p'nx ) 

Jdti  \ 

fdu'  dv,s 


+  (v  +  Vt)\  -  u 


,  fdvf_  dvf_ 
\dx  +  dx 


/  /  /  \®u*  *  / 
nx  +  (unx  +  unx)—  -  u  { 


.  .  5n*  *  f  du'  dw'\  .  ,  ,  \9u 

+  (uriy  +  v  nx)— -  -  u  I  —  +  —  \nz  +  (unz  +  w  nx) 


dy 


dz 


\dy  dx  J  v 
dxdt ,  (2.22) 


where  the  integral  is  understood  to  be  a  surface  integral  over  dQ.  When  added  to 
the  boundary  terms  arising  from  v,  w,  and  p,  the  full  boundary  term  is  given  by 


where  the  last  term  comes  from  the  continuity  equation.  Now,  (2.19)  is  complete. 
The  adjoint  equations  for  q*  are  given  by 


■a/iv  =  r 


with  initial  and  boundary  information 


Qt=r, sen  ~  Jo 
g*(t,x,q*  :  x  G  dfl)  =  0. 


(2.24) 


(2.25) 

(2.26) 


The  adjoint  boundary  conditions,  g*,  depend  on  the  control  set  and  cost  function 
and  will  be  determined  in  later  chapters.  An  important  piece  of  information  to  note 
is  that  the  “initial”  condition  for  the  adjoint  equations  is  specified  at  the  terminal 
time  T.  This  is  because  the  sign  of  the  time  derivative  changes  in  going  from  the 
linearized  equations  to  the  adjoint  equations.  Hence,  characteristics  of  the  adjoint 
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equations  propagate  backwards  in  time,  not  forwards. 

All  of  the  equations  needed  for  this  work  have  now  been  presented  in  their 
continuous  form.  The  next  chapter  will  detail  how  the  continuous  equations  are 
discretized  for  numerical  solution. 


Chapter  3 

Discretization  of  the  Continuous 
Equations* 


3.1  Introduction 

In  Chapter  2  the  continuous  equations  that  are  of  interest  in  this  work  were  de¬ 
rived.  In  order  to  obtain  solutions  to  the  LES  equations  in  the  high  Reynolds  num¬ 
ber  flows  of  engineering  interest,  it  is  necessary  to  solve  them  numerically.  Since 
adjoint-based  optimization  techniques  are  used,  an  important  issue  that  must  be 
addressed  is  how  the  adjoint  equations  are  discretized.  As  mentioned  in  Chapter 
2,  the  approach  taken  iis  to  first  derive  a  continuous  set  of  adjoint  equations  from 
the  continuous  LES  equations,  and  then  discretize  them.  The  alternative  is  to 
formulate  discrete  adjoint  equations  directly  from  the  discrete  LES  operator.  In 
applications  involving  optimization  of  physical  systems,  it  is  unclear  which  formu¬ 
lation  is  superior.  Adjoints  derived  from  the  discrete  equations  typically  achieve 
greater  cost  function  reduction,  however,  some  of  this  reduction  may  be  unphysi¬ 
cal  in  that  the  optimization  takes  advantage  of  peculiarities  in  the  discrete  system 
that  do  not  exist  in  the  continuous  system.  This  type  of  adjoint  is  also  often  more 

’Portions  of  Chapters  3  and  4  have  been  published  as 
Templeton,  J.A.,  Wang,  M.,  &  Moin,  P.  2006  An  efficient  wall  model  for  large-eddy  simulation 
based  on  optimal  control  theory.  Physics  of  Fluids  18(2),  pp.  025101-1  -  025101-13 
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complex  to  derive  and  evaluate. 

In  the  present  application,  it  is  clear  that  for  maximum  accuracy  to  be  achieved, 
the  adjoint  equations  formulated  from  the  discrete  LES  system  should  be  used. 
This  is  because  the  entire  purpose  of  the  optimization  is  to  manipulate  the  dis¬ 
crete  system,  not  to  optimize  an  engineering  cost  function  in  a  physical  system. 
Unfortunately,  this  observation  conflicts  with  the  need  to  have  the  evaluation  of 
the  adjoint  equations  be  as  inexpensive  as  possible.  It  is  therefore  necessary  to 
derive  the  continuous  adjoint  equations  and  consider  only  a  subset  of  the  possible 
methods  used  for  the  time  advancement.  A  further  complication  is  that  the  choice 
of  the  cost  function  can  impact  the  performance  of  various  temporal  discretization 
schemes.  This  effect  will  result  in  certain  cost  functions  proving  more  appropriate 
than  others  for  this  application.  The  effectiveness  of  the  cost  functions  coupled 
with  the  temporal  discretization  will  be  evaluated  in  Chapter  4  based  on  the  accu¬ 
racy  of  the  prediction  of  the  mean  velocity  profiles,  the  cost,  and  the  potential  for 
approximation  (to  further  reduce  the  computational  cost). 


3.2  Discretization  of  the  LES  Equations 

Many  techniques  exist  for  discretizing  the  Navier-Stokes  equation.  For  an  intro¬ 
duction  to  some  of  them,  see  Moin  (2001).  The  specific  method  chosen  for  this 
work  are  low-order  finite  difference  schemes  to  evaluate  the  spatial  derivatives  with 
Runge-Kutta  time  advancement.  These  methods  have  been  selected  because  of 
their  relative  simplicity,  not  for  their  accuracy,  since  the  end  goal  is  to  generate 
wall  models  for  the  coarse  LES  needed  for  industrial  applications.  Higher  order 
numerical  methods,  particularly  spectral  methods,  are  well  suited  to  the  flow  ge¬ 
ometry  that  will  be  used  in  this  work  (see  e.g.  Moser  et  al.  (1999)).  However,  it 
has  been  observed  by  Baggett  et  al.  (2000)  that  the  numerical  techniques  used  in 
a  simulation  can  affect  what  is  required  wall  stresses,  so  these  will  be  eschewed 
in  favor  of  methods  that  will  be  used  in  engineering  applications  to  evaluate  the 
behavior  of  the  wall  model  in  those  situations. 
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Figure  3.1:  Staggered  grid  schematic. 

The  specific  spatial  discretization  used  will  be  a  centered,  second-order  finite 
difference  technique  on  a  staggered  grid.  A  schematic  illustrating  the  system  is 
shown  in  Figure  3.1  in  a  two-dimensional  plane.  The  LES  equations  are  advanced 
in  time  using  a  low-storage,  third-order  Runge-Kutta  technique  (Spalart  et  al, 
1991)  to  advance  the  convective  terms  and  the  Crank-Nicolson  method  for  the 
diffusion  terms.  The  final  aspect  of  the  time  advancement  that  must  be  addressed 
is  the  computation  of  the  pressure.  In  the  present  work,  this  is  handled  using  the 
fractional  step  method  that  allows  for  pressure  to  be  used  as  a  Lagrange  multiplier 
to  ensure  that  the  resulting  velocity  fields  are  divergence  free  (Kim  and  Moin,  1985). 
The  eddy  viscosity  and  wall-boundary  conditions  are  updated  only  at  the  beginning 
of  each  Runge-Kutta  advancement  to  reduce  the  expense  of  the  method  (Le  et  al, 
1997). 


3.3  Discretization  Approaches  for  Adjoint  Equa¬ 
tions 

Given  the  adjoint  equations  formulated  in  a  continuous  setting,  they  must  be  dis¬ 
cretized  so  that  they  may  be  solved  numerically.  As  previously  mentioned,  there 
are  several  objectives  in  analysing  the  discretization  process,  and  one  of  them  is 
amenability  to  approximation.  Making  the  approximations  necessary  to  have  a  sys¬ 
tem  that  can  be  solved  relatively  easily  will  depend  on  the  form  of  the  equations. 
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Two  approximations  will  be  made  in  all  discretization  attempts.  The  first  is  that 
the  optimal  control  problem  will  be  solved  independently  over  each  time  step.  In 
this  sense,  we  are  computing  a  sub-optimal  control  since  a  fully  optimal  control 
would  require  solving  the  adjoint  equations  over  a  sufficiently  large  time  horizon 
that  the  flow  could  converge  to  a  statistically  steady  state  and  statistics  could  be 
measured.  In  addition  to  the  computational  burden  this  entails,  the  memory  re¬ 
quirements  can  also  be  quite  prohibitive  due  to  the  need  to  store  the  time  history 
of  q  over  this  entire  interval  (Bewley  et  al. ,  2001).  The  second  approximation  made 
will  be  to  use  a  single  step  method  for  time  advancement,  as  opposed  to  the  three 
step  Runge-Kutta  scheme  used  to  solve  the  LES  equations. 

The  first  step  of  the  discretization  will  be  to  consider  the  spatial  derivatives. 
In  all  cases,  they  will  be  approximated  in  a  way  consistent  with  the  second-order 
treatment  of  these  derivatives  in  the  Navier-Stokes  equations.  The  justification  for 
such  a  discretization  is  that  the  derivative  operators  acting  on  equation  (2.14)  are 
matrix  operators  for  spatially  discrete  variables  and  the  adjoint  of  this  matrix  is 
simply  its  transpose.  On  any  grid,  this  means  that  the  two  operators  are  of  the 
same  form  (since  they  are  square),  but  on  a  uniform  grid,  the  matrix  is  symmetric 
and  so  applying  the  same  discretization  is  consistent  with  this  formulation. 

The  more  challenging  discretization  to  construct  is  the  temporal  one.  It  would 
be  possible  to  construct  the  exact  adjoint  of  the  discrete,  three-step  system,  but 
using  a  three-step  approach  is  computationally  burdensome.  Using  a  one-step  ap¬ 
proach  immediately  reduces  the  computational  cost  of  each  iteration  by  2/3.  Fur¬ 
ther,  such  an  approach  will  make  it  easier  to  use  this  technique  in  other  codes 
with  different  time-stepping  schemes  since  it  is  not  desirable  to  have  to  re-derive 
the  adjoint  equations  for  each  new  code.  A  further  complication  in  choosing  the 
time- advancement  approach  is  how  it  interacts  with  the  choice  of  the  cost  func¬ 
tion.  Therefore,  we  first  examine  the  impact  this  choice  has  on  the  discrete  adjoint 
equations. 
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3.3.1  Cost  Function  Options  and  the  Resulting  Adjoint  Sys¬ 
tems 

In  this  section,  different  cost  function  formulations  will  be  presented  and  their  effects 
on  the  discrete  adjoint  equations  analyzed.  The  two  main  formulations  of  the  cost 
function  we  will  consider  here  are  distinguished  according  to  whether  they  measure 
quantities  from  the  LES  at  all  times  or  just  the  terminal  time.  When  considering 
the  optimization  over  one  time  step,  it  may  seem  that  this  is  trivial,  but  it  does 
bring  about  changes  in  the  discrete  adjoint  equations  that  must  be  examined  to 
determine  which  will  produce  the  best  results.  The  cost  function  defined  over  all 
time,  denoted  by  Ji,  is: 


Ji  =  (j(q(t,x)),j(q(t,x))) ,  (3.1) 

where  j  is  a  function  that  maps  the  state  to  the  set  of  of  real  numbers.  This 
inner  product  formulation  is  the  same  as  was  used  to  define  the  adjoint  operator  in 
Chapter  2.  This  makes  J\  a  positive  semi-definite  function. 

Using  this  formulation,  boundary  conditions,  initial  conditions,  and  source  term 
of  the  adjoint  system  can  be  set  such  that  the  gradient  may  be  identified  through 
the  solution  of  the  adjoint  equations.  The  gradients  as  functions  of  the  adjoint 
state  will  be  found  via  the  {M'qq',q*)  term,  the  boundary  values,  or  the  terminal 
(t  =  0)  values  of  the  adjoint  state,  depending  on  whether  the  control  inputs  are 
body  forces,  boundary  conditions,  or  initial  conditions,  respectively.  This  means 
that  the  gradient  of  the  cost  function  must  appear  in  (2.19)  through  either  the 
initial  adjoint  conditions  or  the  source  term  since 


W.Ki')  =  <«',/*>• 
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In  order  to  determine  exactly  how  the  cost  function  gradients  should  be  repre¬ 
sented  in  the  adjoint  identity,  it  is  necessary  to  first  take  the  Frechet  derivative  of 
(3.1)  using  the  chain  rule: 


^  =  2 
D<PV 


■  Dj  Dq~ 

3  Dq'Djf 


(3.2) 


recalling  </>  is  an  arbitrary  function  in  £2((Q  x  (0,T])4).  We  note  though  that  since 
Dj/Dq  is  a  linear  operator  acting  over  all  the  values  of  q1, 


(3.3) 


we  can  rewrite  (3.2)  as 


(3.4) 


This  formulation  suggests  that  the  correct  way  to  incorporate  this  information  into 
the  adjoint  identity  is  to  prescribe  the  adjoint  source  term  to  be 
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Du 

Pi 
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Dj  ’ 
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iDpJ 


while  the  initial  conditions  are 


(3.5) 


=  T,  x)  =  0. 


(3.6) 


This  choice  results  in 


v-w  =  T&- 


An  alternative  to  computing  a  cost  function  that  maintains  a  running  track  of 
the  deviation  of  the  LES  from  its  target  is  to  have  one  that  only  measures  this 
deviation  at  the  terminal  time  T.  Mathematically,  this  cost  function  is  constructed 
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as 

J2  =  1 1 1  [j{q{T'S))]2  (3'7) 
where  j  is  the  same  as  in  3.1.  From  a  physical  perspective,  J\  is  more  appropriate 
since  the  LES  quantities  that  are  controlled  are  well  defined  only  after  plane  and 
temporal  averaging.  However,  since  we  will  only  be  considering  the  sub-optimal 
approach  applied  over  one  time  step,  the  number  of  samples  available  in  each 
framework  will  be  identical.  It  will  therefore  remain  to  be  seen  which  approach 
can  produce  the  most  accurate  results  at  the  least  expense. 

This  formulation  will  have  the  gradients  as  functions  of  the  adjoint  state  identi¬ 
fied  in  the  same  way  as  J\.  Also,  while  the  boundary  conditions  of  the  two  systems 
will  be  the  same,  the  gradient  of  3.7  must  be  evaluated  to  determine 

=  /  /  J  2MT,i))^q(T,x)' dx.  (3.8) 


Since  the  term  depending  on  q'  contains  no  integration  in  time,  the  initial  conditions 
will  be  used  to  generate  the  correct  cost  function  gradient.  Thus,  they  are  taken  to 
be: 


Q  |t=T 


‘2 MT))g 
2 MT))& 

2 i(q(T))§i 
.2 MT))%_ 


(3.9) 


Then  the  correct  adjoint  source  term  is 


f,\  =  0. 


(3.10) 


This  choice  of  initial  and  boundary  conditions,  as  well  as  the  source  term,  results 
in  the  initial  conditions  in  (2.19)  becoming 

«\T  =  J  J  j  v!i(T,x)u*(T,x)dx  =  J  J  J  q'(T)2j(q(T))^dx.  (3.11) 

This  result  is,  not  unexpectedly,  similar  to  the  previous  result  for  J\.  However,  the 
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formulation  of  these  problems  will  make  them  more  or  less  amenable  to  different 
temporal  discretization  schemes,  as  will  be  discussed  next. 


Effect  of  Cost  Function  on  Temporal  Discretization 


The  time  integration  and  differentiation  found  in  the  expressions  for  the  adjoint 
equations  and  the  cost  function  gradients  now  need  to  be  addressed  when  dis¬ 
cretizing  the  equations.  The  time  advancement,  since  we  restrict  our  attention  to 
single-step  methods,  can  be  any  combination  of  implicit  or  explicit  methods  applied 
to  each  term.  Quadrature  for  integration  can  be  handled  similarly.  Determining 
a  good  way  to  do  this  is  important  since  nearly  all  of  these  approaches  will  be 
first-order  accurate  in  time  and  so  we  do  not  know  a  priori  which  will  perform 
the  best  in  this  application.  Further,  the  best  method  may  depend  on  how  the 
Navier-Stokes  equations  themselves  are  advanced.  We  will  now  present  the  equa¬ 
tions  using  several  temporal  discretization  schemes.  These  will  be  compared  in  the 
next  section. 

The  first  time-advancement  approach  we  shall  consider  is  a  fully  explicit  scheme. 
Here,  all  the  spatial  derivatives  and  source  terms  are  evaluated  at  the  beginning  of 
the  time  step.  We  denote  variables  at  the  beginning  of  a  time  step  by  n  and  those 
at  the  end  of  the  time  step  by  n  +  1.  Further,  let  At  be  the  time-step  size  for  a 
given  advancement,  and,  where  appropriate,  the  multiplier  2 (3  will  be  used  in  front 
of  it  when  we  wish  to  consider  the  changes  over  only  one  Runge-Kutta  sub-step. 
This  is  consistent  with  the  time  advancement  approach  used  for  the  Navier-Stokes 
equations.  Therefore,  an  explicit  scheme  will  be 


u. 


i,n+ 1 


*  ,  A  J  9utn  ,  9uln 

~  Ui,n  +  At  yUj,n  dx_  +  %n  dx_ 


dp*n  d 

+  +  s; + v,) 


(3.12) 


The  physical  variables  are  aligned  with  the  adjoint  ones  by  noting  that  the  time  of 
T  for  the  physical  variables  corresponds  to  time  0  of  the  adjoint  state  and  n  =  0 
denotes  the  initial  conditions  of  the  adjoint  equations.  In  all  cases  considered,  the 
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adjoint  pressure  will  be  evaluated  using  a  fractional  step  formulation  whereby  an 
intermediate  value  for  u*n+1  is  computed  using  every  term  from  equation  (3.12) 
except  for  the  adjoint  pressure  gradient.  The  adjoint  pressure  is  then  computed 
in  the  same  way  as  the  physical  pressure  by  solving  a  Poisson  equation  with  the 
appropriate  source  term  such  that  subtraction  of  the  pressure  gradient  from  the 
intermediate  adjoint  state  removes  its  divergence.  In  practice,  fully  explicit  schemes 
such  as  this  are  rarely  used  since  they  are  unstable.  However  in  this  application, 
since  we  will  only  use  the  advancement  over  one  time  step,  this  issue  is  not  relevant 
and  so  any  method  can  be  applied. 

Next  an  implicit  time  advancement  scheme  is  considered.  In  this  case,  all  the 
terms  are  evaluated  at  time  level  n  +  1: 

Ui,n+1 
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d<n+l  duj  n+1 


To  reduce  computational  expense,  vt  is  not  evaluated  at  time  level  n+1  which  is 
consistent  with  the  assumption  that  the  control  cannot  significantly  affect  its  value. 
By  incorporating  the  parameter  ip  E  [0, 1],  we  can  generalize  the  above  methods  by 
writing 
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(3.14) 


By  selecting  a  value  for  ip,  the  resulting  expression  can  be  weighted  to  any  degree 
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towards  an  explicit  or  implicit  formulation.  This  can  also  be  done  individually 
term  by  term.  Note  that  the  advancement  for  the  Navier-Stokes  equations  used  is 
either  fully  explicit  or  fully  explicit  in  all  terms  but  wall-normal  diffusion,  which  is 
handled  using  a  semi-implicit  scheme  with  ip  —  1/2. 

What  is  important  in  this  application  is  not  the  specific  time-advancement 
scheme  used,  but  how  it  interacts  with  the  temporally  discretized  cost  function 
and  the  approximations  that  have  already  been  made.  Since  the  cost  functions  are 
reduced  to  only  being  evaluated  at  one  instant  in  time,  there  is  no  choice  as  to 
which  integration  quadrature  should  be  used  to  evaluate  them.  Since  the  state  at 
the  initial  time  is  fixed,  the  quadrature  for  quantities  involving  the  state  must  be 
such  that  the  only  contribution  to  temporal  integrals  come  from  the  next  time  step. 
Similarly,  since  only  one  control  input  will  be  found,  integrals  involving  the  control 
inputs  are  evaluated  such  that  the  only  term  in  them  comes  from  the  current  con¬ 
trol  under  consideration.  When  fully  explicit  advancement  for  the  Navier-Stokes 
equations  are  used,  the  control  is  taken  to  be  at  the  same  time  as  the  initial  state 
and  so  the  quadrature  for  the  state  variables  and  control  variables  are  reversed.  For 
the  semi-implicit  formulation,  this  quantity  must  be  taken  as  the  appropriate  time 
average,  and  the  quadrature  adjusted  accordingly. 

Given  these  definitions  of  quadrature,  it  is  seen  that,  when  considered  over  only 
one  time  step  with  appropriate  constants,  J\  =  However,  the  discrete  equations 
used  to  solve  for  their  gradients  are  different.  The  reason  for  this  is  the  order  At 
difference  that  comes  from  the  time  advancement  of  the  initial  state.  As  At  — >  0, 
this  difference  goes  to  zero  so  that  in  the  limit  the  two  gradients  are  identical.  Thus, 
the  deviation  can  be  seen  as  caused  by  approximating  the  functions  discretely  in 
time,  however  such  deviations  can  be  large  since  in  the  LES  we  are  interested  in, 
large  time  steps  will  be  used.  It  is  also  important  to  note  that  while  we  regulate  the 
indicated  cost  functions,  the  quantities  that  are  actually  of  interest  are  the  long¬ 
time  averaged  flow  statistics.  Thus,  we  must  determine  which  formulation  produces 
the  best  results  for  these  quantities. 
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The  difference  in  the  equations  can  be  immediately  seen  when  considering  the 
explicit  discretization  of  each.  For  J1(  the  explicit  evaluation  is: 
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(3.15) 


Note  that  since  the  initial  conditions  are  all  zero,  the  only  remaining  term  is  the 
source  term.  It  is  evaluated  at  time  T  not  because  of  the  explicit  formulation,  but 
to  be  consistent  with  the  cost  function.  Also,  At  —  T  in  this  case,  which  is  why  it 
is  cancelled  out  in  the  denominator  of  the  source  term.  The  pressure  will  also  be 
absent  here  if  the  source  term  is  divergence  free,  otherwise  it  must  be  included. 

Now,  we  can  compare  (3.15)  with  the  same  equation  formulated  from  J2: 
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In  this  case,  all  the  terms  are  retained  because  the  cost  function  information  enters 
through  the  initial  conditions.  Hence,  it  is  likely  that  (3.16)  will  provide  a  better 
gradient  estimate  than  (3.15)  since  it  incorporates  the  turbulent  features  (through 
terms  containing  u)  and  more  information  about  the  structure  of  the  equations 
(through  the  differential  operators). 


46 


CHAPTER  3.  DISCRETE  EQUATIONS 


A  somewhat  different  pattern  emerges  when 
tion.  Here,  for  the  J\  adjoint  equations  we  have 
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Ui,n+1  ~ Ui,n  4"  ^  ( uj,n+l  +  uj,n+l 

+  If  +  +  (^f  +  ^Ir)  +  *<«■<*  = 


f  ,At(u  dUin+l  duj  n+1 

D±  +Al  j’n+1  dxj  +  j'U+1  dXi 
Dw  '  J 


L  Dp  J 

,  ^Pn+l  ,  JL  /  .  v  )  ( dujn+l  ,  9uj,n+ 

dxj  dxj  t  \  dxj  dxi 


while  for  J2  we  have 

,  a  J..  dutn+ 1  ,  ..  dUj,n+ 1 

Ui,n-hl  — ui,n  +  uj,n+l  o  +  uj,n+l  o 


dp*n+i  d  . 

+  +  + 


£1  +A  l  %n+1  dxj  j'n+1  dXi 

Dw  '  3 


+  +  ^- (^  +  ^)  (^|=±i  +  ^±i)  ) . 

OTj  OXj  \  OXj  OXi  J  J 


These  equations  are  identical.  For  semi- implicit  methods,  the  differences  will  be 
found  in  the  explicit  terms,  in  which  J\  will  have  none  while  J2  will  include  them. 
Therefore  for  single  step  methods,  J2  is  more  robust  with  respect  to  the  numerical 
discretization.  To  illustrate  these  results,  numerical  experiments  have  been  con¬ 
ducted  using  the  different  discretization  formulations  and  cost  functions  in  plane 
channel  flow  (see  Chapter  4  for  more  details).  The  predicted  mean  velocity  profiles 
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are  presented  in  Figure  3.2.  Clearly,  both  formulations  perform  well  when  con¬ 
sidered  with  implicit  temporal  discretization,  while  when  explicit  discretization  is 
used,  the  mean  velocity  profile  is  poorly  predicted. 


Figure  3.2:  Effects  of  adjoint  discretization  and  cost  function  choice  for  ip  =  0  (left) 

and  =  1  (right), - :  ti+  =  2.411og(y+)  +  5.2,, - :  Jj, - :  J2> .  : 

single  Runge-Kutta  sub-step  adjoint  evaluation. 

The  next  question  is  that,  if  the  implicit  formulations  produce  the  exact  same 
equations  for  cost  function  we  know  to  be  identical,  would  that  not  be  the  best 
formula  to  use?  The  answer,  at  least  in  some  cases,  is  no.  This  is  because  the 
formulation  of  the  Navier-Stokes  equations  must  be  taken  into  account.  The  so¬ 
lution  of  the  Navier-Stokes  system  is  computed  using  a  three-step  Runge-Kutta 
advancement  with  each  sub-step  either  being  a  fully  explicit  advance  or  a  semi- 
implicit  advance  with  wall  normal  diffusion  handled  using  ip  =  1/2  and  all  other 
terms  being  fully  explicit.  In  these  two  cases,  it  will  be  straightforward  to  derive 
the  appropriate  adjoint  equations  exactly  from  the  temporally  discretized  state 
equations. 

We  will  first  consider  the  semi-implicit  case  since  that  fits  into  the  above  frame¬ 
work  more  easily.  Here,  the  equation  for  the  Frechet  derivative  of  v*  is 


(3.19) 
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This  will  lead  to  the  adjoint  equations  becoming 


(3.20) 


In  this  case,  very  little  information  from  the  turbulent  flow  is  included  in  this  for¬ 
mulation  since  only  the  initial  conditions  and  vt  depend  on  the  LES  velocity  field. 
Since  Nicoud  et  al.  (2001)  demonstrated  the  importance  of  prescribing  the  wall 
stress  fluctuations,  it  is  important  that  the  turbulent  state  enter  the  adjoint  equa¬ 
tions  so  that  the  resulting  wall  stresses  contain  this  unsteadiness.  Direct  evidence 
of  this  is  also  provided  in  Fig.  3.2  since  the  mean  velocity  profile  computed  using 
this  formulation  is  almost  indentical  to  the  explicit  adjoint  construction  using  J\. 


The  final  formulation,  when  a  fully  explicit  time  advancement  scheme  is  used, 
will  be  demonstrated  to  provide  even  less  information  from  the  turbulent  state. 
Unfortunately,  the  lack  of  derivatives  applied  to  u"+1  means  that  the  previous  for¬ 
mulation  will  not  work  and  we  must  consider  the  fully  discrete  equations.  Before 
doing  this,  we  must  define  the  operator  A,  which  will  be  a  discrete  difference  oper¬ 
ator  in  the  ith  direction: 


^  ®(*^i+i  jXjjXk)  a(xi,Xj,xk) 


(3.21) 


with  each  x  being  an  arbitrary  discrete  coordinate  in  one  direction,  and  hence  the 
need  for  three  of  them  to  specify  a  spatial  location.  The  state  equations  in  this 
formulation  become 


^t,n+l,in  A  iP  =  RHS„ 
un+1,w  -  A up  =  RHS„  +  sign(w)^ 
N(q)=  <  vn+iiW  =  0 

Wn+i,w  ~  A wp  =  RHSn  +  sign(w)!^ 

AjUj,n+l  —  0, 


(3.22) 
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where  “w”  stands  for  quantities  in  cells  adjacent  to  the  walls  and  “in”  stands  for 
quantities  in  the  interior  of  the  domain  (not  adjacent  to  the  walls).  The  terms 
4>u  and  <pw  are  the  streamwise  and  spanwise  control  inputs,  respectively,  which 
enter  (3.21)  in  the  specific  form  they  do  because  the  represent  the  streamwise  and 
spanwise  wall  stresses.  This  occurs  since  only  the  change  in  the  wall  stress  can 
effect  the  state  either  directly  on  the  terms  at  the  wall  (w)  or  through  altering  the 
divergence  and  changing  p  everywhere. 


Taking  the  Prechet  derivative  of  (3.22)  gives  us  the  linear  system 


<n+l,in  -  AiP'  =  0 

<+i,w  “  A„p'  =  sign(w)^ 

<+i,w  =  0 

<+i,w  -  &wP'  =  sign(w)^ 


,n+l 


=  0. 


(3.23) 


In  this  case,  the  inner  product  used  is  a  weighted  summation  that  is  analogous 
to  integration  in  the  continuous  sense.  Since  we  are  on  a  uniform  grid,  the  trans¬ 
pose  operation  to  construct  the  adjoint  equations  can  be  expressed  with  the  same 
notation 


v 


<n+l,in  + 

^n+l,w 

+  A  iP* 

^n+l,w 

=  0 

<+l,w 

-  A  wp* 

—Ajtij 

,71+1  ' 

(3.24) 


Using  the  same  source  term  as  in  the  the  J\  formulation,  we  arrive  at  an  expression 
for  the  gradient  via  the  adjoint  identity 


DJ 

D<t> 


4>  =  2l,J2  <sign(w) 

W 


4>i 

Ay 


Ax  A z, 


(3.25) 
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DJ  .  u\ 

m =  s,gnW  v 


(3.26) 


Note  that  in  the  two  expressions  that  are  exact  for  the  Runge-Kutta  sub-step,  very 
little  information  about  the  physics  or  state  of  the  flow  enters  into  the  formulation. 
It  is  therefore  reasonable  to  expect  that  over  an  entire  Runge-Kutta  step  (or  the 
whole  computation),  these  will  not  perform  as  well  as  the  formulations  starting 
from  the  continuous  equations  since  the  state  and  physics  would  enter  into  the 
adjoint  of  the  full  Runge-Kutta  step. 

While  the  two  cases  considered  above  are  extreme,  it  should  also  be  pointed 
out  that  if  any  single  step  time  advancement  is  used  with  the  parameter  ip  for 
the  Navier-Stokes  equations,  then  the  corresponding  exact  adjoint  equations  over 
one  time  step  will  have  the  same  formulation  but  with  parameter  1  —  ip  on  the 
t  +  At  terms  and  0  on  the  terms  at  time  t.  This  is  because  the  terms  at  time 
t  are  fixed  and  are  therefore  completely  insensitive  to  the  control.  In  this  case, 
it  can  be  seen  that  if  the  Navier-Stokes  equations  are  solved  with  a  fully  implicit 
method,  the  corresponding  adjoint  equations  will  be  fully  explicit,  while  when  the 
Navier-Stokes  equations  are  solved  with  a  fully  explicit  method,  none  of  the  adjoint 
equations  discussed  above  (except  3.24)  are  appropriate.  In  fact,  the  fully  implicit 
adjoint  solution  does  not  correspond  to  any  time  advancement  scheme  over  one 
time  step  for  the  Navier-Stokes  equations  that  does  not  alter  the  terms  at  time  t. 


3.4  Conclusions 

The  analysis  performed  in  this  chapter  can  be  used  to  make  several  choices  with 
respect  to  the  cost  function  and  adjoint  equations  that  will  be  used  in  the  remainder 
of  this  work.  First,  adjoint  equations  based  on  one  Runge-Kutta  sub-step  of  the 
exact  discrete  system  have  been  shown  to  contain  very  little  information  about  the 
turbulent  flow.  The  effect  is  that  they  will  be  ineffective  for  computing  the  fluc¬ 
tuating  component  of  the  wall  stress  required  to  obtain  an  accurate  mean  velocity 
profile.  Therefore,  the  three  step  Runge-Kutta  method  will  be  approximated  by  a 
singe-step  method  over  the  same  interval  instead  of  computing  the  exact  optimal 
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control  for  one  Runge-Kutta  sub-step  and  applying  it  over  all  three  sub-steps. 

Second,  when  using  the  terminal  valued  cost  function  J2,  defined  by  (3.7),  the 
discretization  of  the  adjoint  equations  using  any  single-step  method  retains  a  sig¬ 
nificant  amount  of  information  regarding  the  turbulent  velocity  field.  This  result 
is  in  contrast  to  the  cost  function  J\  defined  using  an  inner  product  (see  (3.1)). 
The  latter  loses  information  about  the  LES  state  as  it  becomes  increasingly  explicit 
through  the  parameter  ip  in  (3.14).  Therefore,  J2  will  be  used  throughout  the  re¬ 
mainder  of  this  work.  The  advantage  of  this  formulation  is  that  a  range  of  time 
integration  methods  can  be  used  with  it. 

The  optimization  procedure  used  in  this  work  is  given  by  the  following  algo¬ 
rithm: 

1.  Compute  the  eddy  viscosity  using  the  dynamic  Smagorinsky  model. 

2.  Forward  solve  the  LES  state  (2.6)  using  an  Euler  step. 

3.  Compute  j  from  (3.7). 

4.  Backward  solve  the  adjoint  equations  (4.13)  to  compute  the  cost  function 
gradient. 

5.  Update  the  control  using  steepest  descent  (4.27). 

6.  Repeat  from  2.  until  the  control  converges. 

7.  Advance  the  LES  state  (2.6)  using  third-order  Runge  Kutta. 

The  next  chapter  will  demonstrate  how  the  computational  cost  of  solving  the  ad¬ 
joint  equations  can  be  reduced  through  a  judicious  choice  of  the  time  advancement 
scheme. 


Chapter  4 

Methods  to  Reduce 
Computational  Expense 

4.1  Introduction 

The  end  goal  of  this  work  is  to  provide  a  general  wall  model  that  can  be  easily  used 
with  an  arbitrary  LES.  In  order  to  achieve  this,  the  three  components  of  the  adjoint 
solver  must  be  relatively  automatic  in  order  to  interface  with  any  given  code  and 
to  be  of  relative  ease  of  use  to  practioners.  These  three  components  are  1)  the 
analytic  derivation,  2)  the  numerical  implementation,  and  3)  the  computational 
resources  required.  The  work  presented  in  discussed  Chapter  3  has  the  advantage 
that,  since  the  adjoint  equations  are  derived  in  their  continuous  form  and  then 
discretized  in  time  and  space,  their  analytic  derivation  is  general  and  need  not  be 
repeated  for  every  code.  The  formulation  also  allows  for  researchers  to  examine 
different  discretization  techniques. 

However,  the  adjoint  solutions  require  an  extensive  programming  effort.  Specif¬ 
ically,  all  the  data  structures  and  routines  to  handle  the  full  set  of  adjoint  equa¬ 
tions  must  be  written.  The  complexity  of  this  code  is  almost  that  of  the  standard 
Navier- Stokes  solver  itself.  Also,  the  adjoint  formulation  depends  directly  on  the 
Navier-Stokes  solver’s  numerical  methods,  primarily  the  spatial  discretization.  It 
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is  impractical  then  to  require  that  a  new  adjoint  solver  be  written  for  every  code 
in  which  a  wall  model  is  desired,  particularly  those  using  more  complex  numerical 
techniques  such  as  curvilinear  coordinates  or  unstructured  meshes.  Preferably,  a 
model  would  be  developed  that  can  be  called  from  a  subroutine  and  used  directly 
with  as  little  modification  as  possible. 

In  addition  to  the  difficulty  involved  with  writing  the  computer  programs  to 
solve  the  adjoint  equations,  their  actual  solution  requires  access  to  significant  com¬ 
putational  resources.  This  is  at  odds  with  the  objective  of  making  LES  easily  af¬ 
fordable.  The  reason  for  this  expense  is  that  the  adjoint  equations  require  roughly 
the  same  computational  effort  to  solve  as  the  Navier-Stokes  equations.  While  typ¬ 
ically  one-step  methods  are  used  to  advance  these  equations,  when  combined  with 
the  iterations  required  for  the  optimization  as  well  as  the  solutions  of  the  Navier- 
Stokes  equations  at  each  iteration  to  provide  the  adjoint  coefficients,  the  cost  can 
become  prohibitively  expensive. 

In  order  to  work  toward  the  goal  of  affordable  LES,  the  issues  of  algorithmic 
complexity  and  computational  expense  must  be  addressed.  This  chapter  examines 
methods  of  mitigating  both  of  these  problems  by  restricting  the  solution  of  the 
adjoint  equations  to  the  near  wall  region.  The  justification  for  this  restriction  is 
that  it  is  the  near-wall  flow  that  is  where  the  SGS  models  have  the  most  difficulty 
and  where  numerical  errors  are  greatest  due  to  the  large  magnitude  of  the  physical 
second  derivative  of  velocity.  Furthermore,  since  the  adjoint  equations  are  solved 
over  only  one  time  step,  the  sensitivity  of  the  control  to  the  outer  flow  is  small. 
These  observations  will  allow  the  entire  optimization  process  to  be  conducted  only 
over  the  near- wall  region,  reducing  the  complexity  and  expense  of  the  wall  model. 

In  this  chapter,  the  specific  problem  under  consideration  will  be  presented,  and 
a  quick  summary  of  the  application  of  the  general  results  found  in  Chapters  2  and  3 
will  be  applied  to  it.  The  problem  is  that  of  flow  in  a  pressure- gradient  driven  plane 
channel,  and  it  is  chosen  because  of  the  large  body  of  knowledge  available  for  this 
flow.  This  makes  it  much  easier  to  determine  what  exactly  the  wall  model  is  doing 
and  to  obtain  a  deeper  understanding  of  the  dynamics  of  the  system.  In  addition, 
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as  has  been  shown  in  Chapter  1,  this  is  a  challenging  test  case  for  standard  wall 
models  and  the  work  can  be  readily  compared  to  the  previous  work  in  control-based 
wall  modeling  by  Nicoud  et  al.  (2001)  and  Baggett  et  al.  (2000). 


4.2  Application  of  the  Adjoint  Problem  to  Chan¬ 
nel  Flow 


In  this  section,  the  specific  LES,  linearized  LES,  and  adjoint  operators  for  plane 
channel  flow  will  be  extended  from  the  results  in  Chapter  2.  Of  important  note  here 
will  be  the  precise  definition  of  the  applied  control,  the  boundary  conditions,  initial 
conditions,  and  actual  equations.  We  start  with  the  continuous  LES  operator, 


Af{q)  = 


duj 

at 


+ 


dujUj 

dij 


+  &  ~  afj  (u  +  **)  ( 


duj 

dij 


+ 


du 

dx 


f) 


duj 

dij 


(4.1) 


with  boundary  conditions 


du 

dyn 


v  =  0 


dw 

V-7—  = 

dyn 


(4.2) 

(4.3) 

(4.4) 


In  this  application,  the  boundary  conditions  in  the  streamwise  and  spanwise  direc¬ 
tions  are  the  control  inputs,  <j).  In  pressure-gradient  driven  plane  channel  flow,  the 
LES  equations  are  written  succinctly  as 
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with  the  source  term  representing  the  imposed  mean  pressure  gradient.  Note  that 
the  flow  variables  have  been  scaled  by  the  friction  velocity,  uT,  for  the  velocity,  the 
channel  half-height,  h,  for  length,  and  pu2T  for  pressure.  This  results  in  the  mean 
pressure  gradient  always  being  —1  since  the  velocity  is  non-dimensionalized  with 
uT. 


4.2.1  Continuous  Adjoint  Operator 

To  construct  the  adjoint  equations,  it  is  first  necessary  to  find  the  linearized  LES 
operator.  As  in  Chapter  2,  this  is  accomplished  by  taking  the  Frechet  derivative  of 

my 


w  =  { 


9u' 

at 

du : 

dxi 


+  §£-£<, +  *>(£  +  g), 


(4.6) 


Recall  we  have  ignored  the  sensitivity  of  ut  to  changes  in  <f>.  In  the  same  way,  the 
initial  and  boundary  conditions  for  the  equations  are  found  to  be: 


(hi! 

ua~  = 

oyn 

(4.7) 

v'  =  0 

(4.8) 

Bin' 

u  =  (f,w. 

oyn 

(4.9) 

The  linearized  Navier-Stokes  equations  are 


(4.10) 


It  is  clear  that  (4.7)  cannot  be  solved  since  the  boundary  conditions  contain  the 
unknown  function  (f>,  motivating  the  need  for  the  adjoint  equations. 
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Using  the  analysis  developed  in  Chapter  2,  the  adjoint  equations  are  found  to 
be 


(4.11) 


By  applying  the  periodic  boundary  conditions  in  the  homogenous  directions,  the 
only  boundary  terms  remaining  from  (2.23)  occur  at  the  walls  (y  =  ±1)  and  the 
temporal  boundaries: 


BT  =  u'ju* 


-{v  +  vt)u* 


(^+P) 

\  9xj  J 


i 

+  (u+  Vtju'j 

-1 


\dxj  dy ) 


+pV  n,, 

(4.12) 


where  the  following  notation  is  used: 


a 


l  _ 
-l 


nix,  1 ,2,  t)dxdzdt 


—  1  ,z,t)  dxdzdt. 


Note  that  the  boundary  conditions  =  v'\w  =  0  have  already  been  used  to  reduce 
the  expression.  With  these  terms  known,  the  adjoint  identity  (2.19)  is  complete. 
The  adjoint  equations  themselves  define  the  solution  of  q*  that  satisfies 


Kl'  =  (4.13) 

with  initial  and  boundary  information  given  by 

9*=r,fen  =  Qo  (4-14) 

g*(t,x,q*  :  x  e  dQ)  =  0.  (4.15) 


How  to  prescribe  the  source  term,  initial  conditions,  and  boundary  conditions  with 
the  given  cost  function  will  be  the  subject  of  the  next  section. 
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4.2.2  Cost  Function  Definition  and  Resulting  Boundary  and 
Initial  Conditions 


Based  on  the  results  in  Chapter  3,  the  cost  function  used  will  be  the  one  defined 
at  the  terminal  time.  In  this  case,  it  is  given  by 

&  =  J'(K(yf  +  s'Jy?)dy,  (4.16) 

where 

KSy)  =  \f  /  =  T)  ~  uitREF)dxdz.  (4.17) 

Before  specifying  the  initial  conditions  implied  by  this  cost  function,  it  is  first 
necessary  to  examine  the  gradient  of  J2: 


D<pv 


dx. 


(4.18) 


Since  the  term  depending  on  cf  contains  no  integration  in  time,  the  initial  conditions 
will  be  used  to  generate  the  correct  formulation.  Thus,  they  are  taken  to  be: 


%  = 


A 

0 

A 

0 


while  the  boundary  conditions  are 


v*  =  0 

oyn 


(4.19) 


(4.20) 

(4.21) 

(4.22) 
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Then,  the  correct  adjoint  equation  to  solve  is 


(4.23) 


This  choice  of  initial  and  boundary  conditions,  as  well  as  the  source  term,  results 
in  the  following  form  for  (2.19) 


DJ1~1 
D<t> 


4>  =  -f  (  (  sign(y)(u*0i  +  ro*03)  dx  dy  dt,  (4.24) 

JO  J  J  t/=il 


where  application  of  (4.7)  and  (4.19)-(4.20)  have  reduced  the  boundary  terms  from 
(4.14)  to  those  appearing  in  (4.24).  The  expressions  for  the  gradient  of  J2  are  then 


PJ2 

D(j)u 

PJ2 

Dcj)w 


=  -sign(t/)u* 
=  — sign(y)u>* 


(4.25) 

(4.26) 


4.3  Computational  Domain 

In  order  to  test  the  models  proposed  in  this  work,  a  pressure-gradient  driven  plane 
channel  flow  at  ReT  =  4000  will  be  considered  in  order  to  compare  with  the  LES 
results  of  Kravchenko  et  al.  (1996).  The  boundary  conditions  used  are  periodic  in 
the  streamwise  (re)  and  spanwise  ( z )  directions,  with  wall  stress  conditions  applied 
to  the  streamwise  and  spanwise  velocities  (respectively  u,  w)  at  the  walls  located 
at  y  =  ±1.  The  total  channel  dimensions  are  27t  x  2  x  47t/3,  after  being  made 
dimensionless  by  the  channel  half- height,  h.  The  penetration  velocity  v  at  these 
walls  will  be  set  to  zero.  In  order  to  reduce  the  effort  required  to  solve  the  adjoint 
equations,  the  mean  pressure  gradient  is  kept  constant  while  the  mass  flow  rate  is 
allowed  to  fluctuate. 

A  centered  second-order  finite  difference  scheme  is  used  on  a  staggered  grid  using 
32  x  33  x  32  points  uniformly  distributed  in  each  direction.  Velocity  components 
are  stored  at  cell  faces  to  which  they  are  normal  while  pressure  and  viscosity  are 
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stored  at  the  cell  centers.  Note  that  this  configuration  implies  that  the  data  needed 
directly  at  the  boundaries  are 

du  ,  dw 

uv  —  v—,  v ,  and  wv  —  i/-— 
dy  dy 

which  will  be  provided  by  the  control. 

To  advance  the  equations  in  time,  a  third-order  low-storage  Runge-Kutta  scheme 
is  used.  At  each  Runge-Kutta  sub-step,  the  momentum  equations  are  advanced  ex¬ 
plicitly  except  for  the  wall-normal  diffusion  terms,  which  are  handled  using  the 
Crank-Nicolson  technique.  In  wall-resolved  simulations,  this  is  done  to  avoid  the 
wall-normal  CFL  constraint.  Due  to  the  coarse  grids  used  in  this  work,  this  con¬ 
straint  is  not  as  significant  but  this  method  is  retained  since  it  does  not  substantially 
increase  the  computational  effort.  The  time  step  is  fixed  at  A tuT/h  =  0.0015,  which 
produces  a  maximum  CFL  value  of  approximately  0.3.  Given  the  time  advance¬ 
ment  method  used,  this  is  well  below  the  bounds  required  for  stability,  but  it  was 
found  to  be  necessary  to  have  the  solution  independent  of  At. 

4.3.1  Optimization  Technique 

In  the  optimization  routine,  it  is  necessary  to  iterate  on  the  adjoint  and  LES  equa¬ 
tions  in  order  to  compute  the  cost  function  gradients.  As  previously  discussed  for 
the  adjoint  equations,  a  one-step  method  in  time  will  be  used  to  advance  both  sets 
of  equations  within  the  optimization  routine  in  order  to  reduce  the  computational 
expense  of  the  method.  The  LES  equations  are  advanced  for  one  complete  time- 
step  using  an  explicit  Euler  scheme,  thus  obtaining  the  state  variables  required  in 
(4.11). 

If  an  implicit  solution  is  used  to  solve  the  adjoint  equations,  it  must  have  a 
good  initial  guess  to  avoid  many  iterations.  In  this  work,  it  is  taken  to  be  the 
right-hand  side  vector  of  (3.18),  which  is  an  0(1)  approximation  to  the  solution. 
Note  that  Nicoud  et  al.  (2001)  took  the  initial  guess  to  be  the  adjoint  solution  from 
the  previous  iteration,  which  resulted  in  slower  convergence  of  the  optimization 
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procedure. 

With  the  gradients  determined  from  the  adjoint  state,  the  control  is  updated 
using  a  steepest  descent  method.  Let  the  index  k  denote  the  sub-iterations  per¬ 
formed  in  the  optimization  routine,  and  so  the  update  equation  for  the  controls  (j> 
is 

=  (4-27) 


where  ^  is  a  descent  parameter  set  a  priori  to  be  5  x  105  to  match  the  value  used 
by  Nicoud  et  al.  (2001).  Theoretically,  this  method  is  not  guaranteed  to  converge, 
however  in  this  application  this  has  not  been  observed  to  be  an  issue.  Again,  to 
match  Nicoud  et  al,  convergence  is  measured  in  the  Z^-norm  of  the  change  in 
control  between  sub-steps: 


uk+i  -  <frkr 
ft 


<c 


(4.28) 


with  e  =  2  x  10~5.  In  this  case,  4>s  is  of  order  unity  because  the  mean  wall  stress  is 
one  when  uT  is  used  as  the  velocity  scale.  This  approach  is  recommended  for  cases 
in  which  the  scales  of  the  cost  function  are  unknown  (Dennis,  1983). 

In  order  to  regularize  the  cost  function,  a  penalty  term  on  the  fluctuations  of 
the  controls  is  included.  This  is  done  by  defining  the  cost  function  in  the  following 
manner: 


J  —  Jt  +  Jp,i  +  Jp,3,  Jp,i  =  -jrji  J  j  J(ft- ft)  dxdz.  (4.29) 


This  is  necessary  to  prevent  the  control  from  using  very  large  fluctuations  to  reduce 
the  cost  function,  which  can  destabilize  the  simulation.  In  the  present  work,  a  = 
10-6,  corresponding  to  the  value  used  by  Nicoud  et  al.  (2001).  In  general,  a  should 
be  set  to  the  smallest  value  that  allows  the  simulation  to  be  performed  stably. 
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4.4  Techniques  to  Reduce  Computational  Expense 

4.4.1  Definition  of  a  Near- Wall  Cost  Function  and  its  De¬ 
pendence  on  Pressure  Variables 

An  important  goal  of  this  work  is  to  reduce  the  amount  of  computational  effort  that 
the  wall  modeling  routine  must  exert  away  from  the  wall.  Doing  so  has  two  primary 
advantages.  The  first  is  that  the  operation  count  and  data  storage  requirements  of 
the  model  are  significantly  reduced  if  only  the  data  near  the  wall  must  be  stored 
and  manipulated.  Secondly,  by  defining  the  cost  function  near  the  wall,  it  may  be 
possible  to  use  TBLE  similar  to  Wang  and  Moin  (2002)  to  provide  a  target  for  the 
controller.  This  subsection  illustrates  how  such  a  construction  can  be  accomplished. 

The  first  step  in  restricting  the  equations  solved  in  the  wall  model  to  the  near¬ 
wall  region  is  to  define  a  cost  function  that  only  measures  the  flow  there.  Take  the 
wall  to  be  located  at  y =  0  and  define  y  —  ym  to  be  the  upper  edge  of  the  near- wall 
region.  Then  the  cost  function  is  written  as 

ry-m 

h=  \  (Sl  +  52w)dy.  (4.30) 

Jo 

Unfortunately,  the  smallest  value  of  ym  that  will  produce  a  good  solution  is  un¬ 
known,  and  this  value  could  vary  with  the  Reynolds  number,  numerical  techniques, 
and  grid  resolution  used.  However,  a  reasonable  value  can  be  obtained  by  examining 
the  flow  obtained  using  standard  wall  models.  In  these  predictions,  an  unphysical 
transition  occurs  in  the  mean  profile  from  the  near-wall  flow,  in  which  the  slope  of 
the  logarithmic  profile  is  under  predicted,  to  the  outer  flow,  which  more  accurately 
represents  this  slope.  It  is  reasonable  to  assume  that  it  is  in  this  region  where  SGS 
and  numerical  errors  are  dominant,  and  so  this  is  where  the  cost  function  should  be 
measured.  In  the  present  LES  ( ReT  =  4000),  this  transition  occurs  after  the  third 
wall-normal  velocity  grid  point  at  y  =  0.15/ h  or  y+  =  605,  so  m  is  chosen  to  match 
this  value.  Figure  4.1  shows  the  mean  velocity  profiles  obtained  when  choosing  m 
corresponding  to  the  second  and  third  streamwise  velocity  nodes,  with  figure  4.2 
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showing  the  RMS  velocity  fluctuations.  The  latter  case  shows  significantly  better 
agreement  than  the  former  with  the  results  obtained  when  the  cost  function  is  de¬ 
fined  over  the  whole  domain,  providing  evidence  that  the  supposition  about  the 
transition  region  is  correct.  This  same  region  can  be  determined  in  other  flows  by 
performing  a  coarse  LES  using  a  standard  wall  model  and  noting  where  the  tran¬ 
sition  from  near-wall  to  outer  flow  takes  place,  although  the  results  of  this  work 
suggest  that  the  number  of  overlap  grid  points  is  independent  of  grid  resolution 
and  Reynolds  number. 


Figure  4.1:  Mean  velocity  profiles  at  ReT  =  4000, - :  u+  =  2.411og(y+)  +  5.2, 

- :  full  channel  cost  function, - :  y*  =  605  (3  points), .  :  y+  =  363  (2 

points). 


Now  that  a  cost  function  has  been  defined  near  the  wall,  it  is  important  to 
determine  its  sensitivity  to  the  flow  variables.  The  first  variable  that  is  considered 
will  be  the  pressure.  Since  a  fractional  step  method  is  being  used  to  advance  the 
state  equations,  the  state  can  be  written  as 


Ui  =  u\-  At 


(4.31) 
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Figure  4.2:  RMS  velocity  fluctuations  at  ReT  =  4000, - :  Kravchenko  et  al. 

(1996), - :  full  channel  cost  function, - :  y+  =  605  (3  points),  .  : 

y+  =  363  (2  points). 


with  being  the  intermediate,  non-divergence  free  state  arising  from  advancing  the 
momentum  equations  without  the  most  recent  pressure.  Inserting  this  into  (4.17) 
results  in 

<My)  =  J /KlEs(*,T)  -  A t~-  -  uiiREF(x,T))  dxdz 

=  J  J  (ulES(X,  T)  -  «i,REF(^,  T))  dxdz-  At  J  J  dx  dz. 

Since  the  second  integral  is  over  the  homogenous  directions  and  involves  derivatives 
in  them,  it  is  identically  zero.  Hence,  the  final  evaluation  of  the  pressure  does  not 
have  an  effect  on  the  cost  function  value. 

In  order  to  take  advantage  of  pressure  not  entering  the  cost  function  value,  the 
only  LES  state  data  that  can  be  used  is  that  from  the  old  time.  If  data  at  the 
new  time  is  required,  it  is  coupled  with  the  adjoint  operator  and  can  affect  the 
value  and  distribution  of  the  adjoint  solution.  This  implies  that  the  best  choice  for 
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the  discretization  of  the  adjoint  equations  is  a  fully  implicit  method.  Using  this 
technique,  all  the  LES  state  information  appearing  as  coefficients  in  the  adjoint 
equations  is  at  the  old  time.  This  means  that  a  pressure  solution  for  the  LES 
equations  is  not  required  to  perform  the  optimization.  Then,  only  one  iteration  on 
the  implicit  equations  will  be  performed,  which  is  equivalent  to  the  approach  taken 
by  Nicoud  et  al.  (2001). 


4.4.2  Reduction  in  Expense  of  Navier-Stokes  and  Adjoint 
Solutions 

If  a  one  step,  explicit  advance  is  used  for  the  LES  state  in  the  adjoint  equations,  the 
effects  of  the  control  will  never  be  applied  to  any  point  other  than  the  one  nearest 
the  wall.  However,  if  the  wall-normal  diffusion  is  handled  implicitly,  then  the  effects 
of  the  wall  stresses  will  be  applied  everywhere  in  the  flow.  An  important  property 
that  can  be  leveraged  here  is  the  linearity  of  this  term.  The  time  advancement  can 
be  written  succinctly  as  (with  no  summation  implied  over  the  index  i ) 

Wiui  =  fi  +  bi,  (4.32) 

with  W  being  a  matrix  representing  the  implicit  wall-normal  diffusion  and  time 
advancement,  and  /  and  b  being  the  right-hand  side  vectors  arising  from  the  discrete 
LES  equations  and  boundary  conditions,  respectively.  Since  /  will  not  change 
during  the  course  of  the  iterations,  the  term 

u{  =  W-'fi 


need  only  be  solved  once  per  optimization  procedure.  Since  /  has  non-zero  entries 
everywhere,  this  equation  must  be  solved  over  the  entire  domain.  The  total  velocity 
field  at  the  new  time  is  then  written  as 


Ui  =  u{  +  W{  1bl. 


(4.33) 
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Note  that  since  the  cost  function  only  measures  the  deviation  in  the  streamwise 
and  spanwise  velocity  profiles  from  their  mean  values,  (4.33)  need  not  be  solved  for 
the  wall-normal  velocity  component  in  fully  developed  channel  flow. 

Approximations  can  also  be  made  in  the  solution  of  the  adjoint  equations.  Since 
a  fully  implicit  method  is  used  to  advance  these  equations,  the  only  term  that  will 
appear  on  the  right-hand  side,  and  hence  in  the  initial  guess  for  the  adjoint  state, 
will  be  the  adjoint  initial  conditions.  This  state  varies  only  in  the  wall-normal 
direction  because  it  is  the  plane-wise  averaged  deviation  of  the  LES  and  target 
velocities,  so  many  of  the  terms  that  appear  in  (4.11)  are  zero.  When  only  the 
gradients  in  the  wall-normal  direction  are  retained,  along  with  removing  v*  and  p* 
since  they  are  zero,  the  adjoint  operator  becomes 

Furthermore,  since  the  initial  conditions  do  not  vary  with  x  or  z,  the  gradients  need 
only  be  computed  once  for  each  wall,  which  is  the  advantage  when  using  this  form 
of  the  adjoint  “convective”  terms.  If  the  other  formulation  is  used,  gradients  of  the 
local  turbulent  velocities  must  be  obtained,  making  the  method  more  expensive 
(see  Section  4.2.1)  but  more  generally  applicable. 

4.4.3  Near- Wall  Approximation  to  Implicit  Equations 

In  the  optimization  process,  the  LES  state  only  contributes  to  the  cost  function 
over  the  matching  region  near  the  wall.  Similarly,  the  adjoint  solution  only  affects 
the  control  inputs,  and  hence  the  LES,  through  its  value  at  the  wall.  Therefore, 
the  solution  of  the  equations  throughout  the  outer  domain  contributes  relatively 
little  to  the  overall  performance  of  the  LES.  The  difficulty  in  eliminating  this  outer 
region  is  that  it  does  affect  the  near- wall  region  through  the  implicit  equations  that 


66  CHAPTER  4.  METHODS  TO  REDUCE  COMPUTATIONAL  EXPENSE 


are  solved  to  advance  the  LES  and  adjoint  states.  These  are  the  adjoint  pressure 
equation  and  implicitly  discretized  wall-normal  diffusion.  The  wall-normal  diffusion 
can  be  handled  directly. 

The  structure  of  b  in  (4.33)  can  be  used  to  reduce  the  computational  effort  in 
its  solution.  This  is  done  by  noting  that  the  only  non-zero  entries  of  b  will  be 
those  at  the  wall  containing  the  applied  control  inputs.  The  effect  of  these  terms 
far  away  from  the  wall  over  one  time  step  will  be  quite  small.  Therefore,  W  will 
be  approximated  by  Wnw  which  is  identical  to  W  near  the  wall  but  contains  only 
zeros  away  from  it,  with  a  Neumann  boundary  condition  used  to  approximate  the 
reduction  in  sensitivity  of  the  outer  flow  to  the  wall  stress.  This  condition  will  be 
placed  on  point  m  +  1  so  as  to  mitigate  the  effect  of  this  approximation  on  the 
region  that  defines  the  cost  function.  The  total  velocity  field  at  the  new  time  is 
then  written  as 

u<  =  u{  +  wr!wbt.  (4.35) 

Using  this  technique,  only  the  near- wall  wall-normal  diffusion  needs  to  be  computed 
for  the  LES  system  during  the  optimization  iterations,  significantly  reducing  the 
computational  expense. 

The  final  equation  that  needs  to  be  addressed  is  for  the  pressure  in  the  adjoint 
solution.  Again,  the  only  non-zero  terms  will  be  near  the  wall,  specifically  at  the 
first  m  +  1  points  in  the  wall-normal  direction,  but  the  pressure  can  affect  the  field 
in  the  whole  domain.  In  order  to  advance  these  equations,  the  pressure  is  modified 
in  a  similar  manner  to  the  wall-normal  diffusion  in  the  LES  equations.  A  Neumann 
condition  is  used  to  close  the  adjoint  pressure  solution  at  m+  1.  This  is  reasonable 
since  the  only  pressure  values  that  will  affect  the  controls  will  be  those  at  the  wall 
if  only  one  iteration  is  performed.  Therefore,  as  long  as  the  value  near  the  wall  is 
a  good  approximation  of  the  correct  adjoint  pressure,  the  adjoint  equations  need 
only  be  solved  near  the  wall  as  well.  In  Figures  4.3  and  4.4,  the  adjoint  equations 
are  solved  only  at  the  first  m  +  1  points  and  then  the  pressure  is  computed  as  has 
just  been  described.  The  results  demonstrate  that  the  method  is  still  capable  of 
accurately  predicting  the  mean  velocity  profile. 
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Figure  4.3:  Mean  velocity  profiles, - :  u+  —  2.41  log(y+)  +  5.2, - :  original 

formulation, - :  reduced  cost  formulation. 


Figure  4.4:  RMS  velocity  fluctuations,  -  :  Kravchenko  et  al.  (1996), - : 

original  formulation, - :  reduced  cost  formulation. 
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4.5  Conclusions 

A  method  has  been  presented  that-  significantly  reduces  the  computational  effort 
required  to  implement  a  wall  model  based  on  sub-optimal  control  theory.  Without 
this  method,  the  wall  model  is  approximately  13  times  as  expensive  as  the  LES 
solution,  while  with  the  method,  the  wall  model  is  between  1  and  2  times  the  cost 
of  the  LES.  A  further  advantage  of  this  approach  is  that  only  data  near  the  wall  is 
needed  to  define  the  cost  function,  and  the  optimization  procedure  must  only  solve 
the  LES  and  adjoint  equations  in  this  region.  Thus,  applying  this  technique  to  flows 
in  complex  geometries  is  more  feasible  since  near- wall  approximations  can  be  used 
in  defining  a  predictive  target  profile  and  constructing  the  numerical  techniques 
used  in  the  optimization  routine. 

In  order  to  obtain  such  advantages,  the  structures  of  the  cost  function  and  ad¬ 
joint  equation  discretization  are  exploited.  This  paradigm  could  have  utility  in  other 
applications  where  limited  computer  resources  dictate  that  significant  increases  in 
the  efficiency  of  the  optimization  procedure  is  more  important  than  a  very  accurate 
gradient  evaluation.  In  such  cases,  the  most  general  formulation  of  the  continuous 
adjoint  equations  should  be  found,  which  be  discretized  and  approximated  to  max¬ 
imize  accuracy  while  minimizing  expense.  Therefore  computational  efficiency  can 
be  taken  into  consideration  in  the  construction  of  the  optimization  procedure. 


Chapter  5 

Coupling  LES  and  RANS  via  the 
Sub-Optimal  Control 
Formulation* 


5.1  Introduction 

In  this  chapter  we  will  discuss  how  RANS  equations  can  be  utilized  in  the  optimal 
control  framework  to  provide  the  target  profile  needed  by  the  controller.  There 
are  many  ways  in  which  RANS  can  be  incorporated  with  varying  models,  physics, 
and,  perhaps  most  importantly,  coupling  to  the  rest  of  the  system.  We  will  first 
consider  a  general  RANS  model  arbitrarily  coupled  to  the  LES  state  and  the  con¬ 
trol  inputs.  To  make  the  presentation  more  succinct,  only  the  terminal  time  cost 
function  without  penalties  will  be  considered,  but  the  analysis  could  be  applied  to 
any  formulation.  Penalty  terms  will  be  examined  separately  later  in  the  chapter. 
Also,  the  term  RANS  will  be  used  loosely  to  describe  any  system  that  can  capture 
an  “average” -type  behavior  of  the  near-wall  physics,  which  will  ensure  the  spa¬ 
tially  and  temporally  averaged  LES  solution  accurately  captures  the  mean  velocity 

’Chapter  5  is  currently  being  prepared  for  journal  submission  as 
Templeton,  J.A.,  Wang,  M.,  &  Moin,  P.  A  predictive  wall  model  for  large-eddy  simulation  based 
on  optimal  control  techniques. 
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profile. 

As  a  matter  of  practicality,  it  is  best  to  use  the  cheapest  RANS  system  that 
meets  the  above  requirements.  However,  an  important  constraint  is  that  the  cost  to 
solve  the  RANS  system  can  only  scale  weakly  with  the  Reynolds  number.  This  is 
achieved  by  having  he  RANS  system  resolve  the  strong  wall-normal  mean  velocity 
gradients  while  retaining  the  same  coarseness  in  the  wall-parallel  directions  as  the 
LES.  However,  in  the  derivation  of  the  adjoint  equations  the  RANS  system  will  first 
be  treated  a  general  continuous  operator.  The  specific  RANS  system  used  will  be 
considered  in  the  discretization  process.  Examples  of  several  appropriate  near-wall 
RANS  equations  of  varying  levels  of  complexity  are  presented  in  Wang  and  Moin 
(2002). 

5.2  LES  Control  Algorithm  using  RANS  Targets 

A  control  algorithm  will  be  devised  to  compute  the  wall  stresses  that  are  required 
as  inputs  to  the  LES.  However,  in  order  to  apply  a  control  algorithm,  it  is  necessary 
to  have  a  quantity  to  regulate.  Many  such  flow  quantities  are  possible,  including 
both  volume  integrals  and  wall-normal  fluxes  of  flow  properties  such  as  momentum, 
stress,  energy,  and  vorticity.  However,  since  a  predictive  method  is  required,  it  will 
be  necessary  to  define  a  cost  function  that  can  use  data  obtained  from  RANS-based 
computations.  This  is  because  RANS  is  the  only  technique  that,  at  present,  can 
obtain  reasonably  accurate  information  at  a  modest  cost  in  most  attached  near¬ 
wall  flows.  The  primary  cost  function  considered  will  then  be  based  on  the  plane 
averaged  deviation  of  the  LES  profile  from  the  RANS  profile,  although  other  cost 
functions  will  be  considered.  An  additional  reason  for  using  such  a  cost  function  is 
that  previous  calculations  by  Nicoud  et  al.  (2001)  and  those  presented  in  Chapter 
4  demonstrated  that  this  definition,  at  least  for  the  case  of  a  fixed  target  velocity, 
was  able  to  yield  accurate  mean  velocity  profiles. 

The  RANS  model  will  be  of  the  form  used  by  Wang  and  Moin  (2002),  since 
they  demonstrated  the  effectiveness  of  near- wall  RANS  modeling  in  supplying  wall 
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stresses  for  a  moderately  resolved  LES  in  complex  geometry.  The  cost  function  is 
then  defined  to  be 

J(4>)=  (fiu{y)2 +  5w(yf)dy,  (5.1) 

Jo 

which  is  a  function  of  the  control  input  <p.  In  the  case  of  channel  flow,  (f>  E  (C2(dCl  — > 
R))2  since  each  element  is  a  function  mapping  points  on  the  walls  to  R.  There  are 
two  such  functions  corresponding  to  each  of  t™2  and  t™2  on  both  walls.  For  clarity, 
the  notation  <pu  =  t™2  and  <f)w  =  t™2  will  be  used. 

With  X  and  Z  used  to  designate  the  extent  of  the  domain  in  the  streamwise 
and  spanwise  directions,  the  specific  quantities  appearing  in  the  first  term  of  (5.1) 
are  defined  as 

Sm  (y)  =  -j  J  Jj,uithEs{x,y,  z,  T)  -  in, rans(x,  y,  z,  T))dx  dz,  (5.2) 

where  the  subscripts  l  and  r  denote  state  variables  from  the  LES  and  RANS  equa¬ 
tions,  respectively.  Note  that  both  are  functions  of  <f>,  with  u^les  satisfying  the 
constraint  that  J\f(q)  =  0.  While  the  RANS  velocity  can  be  taken  as  a  fixed  target 
from  an  independent  calculation,  it  is  more  useful  to  have  it  coupled  with  the  LES 
state  and  control  inputs. 


Now  that  the  cost  function  is  defined,  its  derivative  must  be  calculated  to  apply 
the  control  algorithm.  Taking  the  Frechet  derivative  of  J  with  respect  to  (j)  results 
in 


DJ:i  =  2J’" 


D4> 


D(f> 


D<f> 


dy, 


(5.3) 


where  <p  is  an  arbitrary  test  function.  To  evaluate  this  expression,  the  Frechet 
differential  of  5Ui(y )  must  also  be  computed: 


D$ui(y)  7  __  1  f  f  / DUj, LES 

D<f>  A  JxJz\  Dcf> 


■P^t.RANS^ 

D4>  J 


(f)  dx  dz. 


(5.4) 


The  second  term  in  the  integral  can  be  evaluated  using  the  chain  rule  (which  applies 
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for  Frechet  differentiation,  see  Luenberger  (1969)): 

-P^j.RANS  7  _  ( 'PUj.RANS  Z^i.RANS  Dq\  t 

D(j)  *  ~  V  V(f>  Vq  D(j> )  ® 


(5.5) 


where  Vuit rans/T^  represents  the  partial  Frechet  derivative  with  respect  to  (j). 
For  simple  models,  these  derivatives  can  be  computed  analytically  at  virtually  no 
additional  expense.  However,  for  more  complicated  models,  it  will  be  necessary 
to  define  an  adjoint  problem  for  them.  This  procedure  will  now  be  shown  for  a 
general  RANS  model.  The  results  for  specific  model  formulations  will  be  presented 
as  needed. 

Once  the  RANS  sensitivities  have  been  determined,  the  adjoint  source  term  for 
the  LES  adjoint  equations  needs  to  be  modified  to  incorporate  this  information, 
specifically  the  last  term  in  (5.5).  In  order  to  determine  what  the  source  term  for 
the  LES  adjoint  equations  should  be,  it  is  necessary  to  examine  the  structure  of 
the  RANS  sensitivity  operator.  In  the  cases  considered  here,  it  takes  the  form  of  a 
linear  functional  acting  on  q'(x ): 


VUi  rans  / 

— tS - 0 

Vq 


P^i.RANS 

Vqj 


dal. 


(5.6) 


The  only  way  to  have  this  integral  evaluated  within  the  inner  product  formulation 
is  to  have  the  inner  product  of  f*  and  q'  yield  the  expression  in  (5.6).  However, 
the  expression  found  in  (5.3)  and  (5.4)  must  still  be  recovered.  To  this  end,  the 
RANS  gradient  must  be  weighted  by  the  corresponding  term  that  appears  in  the 
cost  function,  which  is  accomplished  by  switching  the  order  of  integration: 
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(f')g'(x')J  da?  )  dydxdz  = 
^t.RANS^)  . 


(a?)  ]  dx  }  dx1. 


(5.7) 


Therefore,  the  adjoint  source  term  will  be  decomposed  as  follows 


/;  =  /*+/;. 


(5.8) 


where  /*  is  the  adjoint  source  term  used  in  Chapter  4  while  /*.  is  given  by 


dx?. 


(5.9) 


The  LES  adjoint  equations  are  then  solved  as  before  to  yield  the  cost  function 
gradient  in  combination  with  the  solution  of  the  RANS  adjoint  equations  for  the 
direct  sensitivity  of  the  RANS  state  to  the  control.  In  this  case,  however,  the  LES 
adjoint  source  term  will  be  a  function  of  all  the  spatial  variables  and  not  just  the 
wall-normal  direction. 


5.2.1  RANS  Sensitivities 

Since  the  RANS  profile  potentially  changes  in  response  to  the  LES  state  and  control 
inputs,  the  sensitivities  of  the  RANS  velocity  profiles  to  these  variables  must  be 
evaluated.  These  sensitivities  are  required  in  order  to  compute  the  necessary  gradi¬ 
ents  that  arise  in  (5.4).  Within  this  framework,  it  is  possible  to  use  many  different 
near-wall  RANS  models,  ranging  from  the  simplest  stress  balance  model  to  a  full 
set  of  time- varying  TBL  equations.  In  addition,  it  is  possible  to  prescribe  a  variety 
of  boundary  conditions  and  source  terms  from  the  LES.  These  include  velocities, 
velocity  gradients,  pressure  gradients,  and  energy  or  vorticity  fluxes  at  ym,  and 
velocities  and  stresses  at  the  wall.  Each  wall  can  then  be  handled  independently. 

In  order  to  compute  the  RANS  sensitivities,  it  is  necessary  to  first  define  a 


74 


CHAPTER  5.  OPTIMAL  CONTROL  FORMULATION  WITH  RANS 


RANS  state,  qr.  Given  the  cost  function  used  in  this  work,  it  is  necessary  that  both 
the  RANS  streamwise  and  spanwise  velocities,  urans  and  wrans>  be  elements  of 
qr.  In  general,  qr  must  include  all  the  states  that  appear  in  the  cost  function  and 
auxiliary  states  which  influence  these  quantities.  For  example,  other  elements  of  qr 
could  be  the  wall-normal  velocity  and  eddy  viscosity  model  variables.  Given  these 
states,  the  RANS  system  can  be  written  generally  as 


Kq(Qr)  =  fr{q,  <P) 


subject  to  boundary  conditions  at  the  wall 


(5.10) 


9qi.9r  |o)  /r,bcw(?)  $)» 


(5.11) 


and  at  ym: 

hq(qr\m)  ~  fr,bcm{Qi4>')-  (5.12) 

Here  the  dependence  of  the  operator,  source  term  and  boundary  conditions  on  q 
is  meant  to  couple  the  RANS  solution  to  the  LES  solution  in  a  potentially  general 
manner.  The  RANS  dependence  on  the  LES  quantities  will  normally  be  only  at 
positions  above  the  matching  layer,  as  they  are  not  considered  to  be  physically 
correct  near  the  wall. 

Now,  taking  the  Frechet  partial  derivative  of  (5.10)-(5.12)  with  respect  to  q  and 
4>  yields  two  equations,  one  for  each  sensitivity: 


Vnq{gr)  VI Zg(qr)  i  _  Vf 

Vq  q  Vqr  Qt  Vqq 
T>TZq(qr)  ~  T>TZ(qr)  $  T)fr~ 

~^Dr^^rq’  =  ^' 


(5.13) 

(5.14) 
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where 


L>qr  ~ 

Vqq 

vy 


Note  that  <fi  and  q  refer  to  perturbations  in  the  LES  state  and  control  input,  re¬ 
spectively,  and  are  taken  to  be  independent. 

With  the  following  definitions 


'RRg(qr) 

Vqr  ’ 

(vfr  _  vng(qr)\  . 

\Vg  Vg  )9, 


(5.15) 

(5.16) 


(5.13)  can  be  written  as  the  linearized  RANS  system 


KA  =  S'Al,  Qr)-  (5.17) 

The  boundary  conditions  for  this  equation  (using  similar  notation  for  the  subscripts) 
are 


4^r|o  =  /;>bcw(g,  qr,cp),  (5.18) 

KaM m  =  /r,bc(9>  Qr),  (5-19) 


where 


/  __  '^9q(.Qr\m) 

3q,qr  - 


f r,bcw(q 


Vqrlm 

_  f  L>fr,  bcw 

Vq 


,9r)  -  ^ 


^ffg(ffrlm) 

Vq 


(5.20) 


(5.21) 


76 


CHAPTER  5.  OPTIMAL  CONTROL  FORMULATION  WITH  RANS 


and 


fr, bc(?)  Qr) 


T>hq(q r|m) 

Hqr  \m 

( Vfr, be  _  Hhq(qr\m)\  _ 
V  Vq  Vq  )Q- 


(5.22) 

(5.23) 


A  similar  system  can  be  derived  for  the  sensitivities  to  <j>. 

Integration  by  parts  now  yields  the  adjoint  operator  R*q  qr  of  (5.17) .  The  adjoint 
identity  then  becomes: 


{I'M,  1r),  «;>  =  «,  /;)  +  BTr.  (5.24) 

The  RANS  adjoint  source  term,  /*,  is  determined  such  that  q'r  can  be  identified. 
Similarly,  the  boundary  conditions  on  q*  must  be  chosen  such  that  all  terms  not 
directly  proportional  to  q'r\m  are  zero.  It  can  be  seen  that  this  will  result  in  a 
well  posed  system  by  considering  that  only  terms  containing  a  derivative  of  or¬ 
der  two  or  greater  in  the  wall-normal  direction  will  require  an  adjoint  boundary 
condition.  This  arises  from  the  integration  by  parts.  Any  terms  containing  first 
order  derivatives  will  have  boundary  terms  proportional  to  elements  of  q'r  and  q*. 
A  derivative  of  order  n  on  an  element  of  q'T  will  result  in  one  term  proportional  to 
q'r  and  n  —  1  terms  proportional  to  derivatives  of  it.  Thus,  all  the  boundary  terms 
involving  derivatives  of  q'r  can  be  removed  by  setting  the  appropriate  derivatives  of 
q*  to  zero.  This  is  also  the  correct  number  of  boundary  conditions  required  for  well 
posedness.  Hence,  equation  (5.24)  allows  for  the  identification,  in  the  weak  sense, 
of  ^  which  can  be  used  to  compute  the  cost  function  gradients. 

5.2.2  RANS  Sensitivities  for  a  Simplified  System 

In  order  to  provide  an  illustrative  example,  as  well  as  to  derive  the  equations  used 
in  this  work,  the  RANS  sensitivities  of  the  previous  section  will  be  applied  to  a 
simple  near-wall  model.  The  model  under  consideration  is  the  simplest  one  used 
by  Wang  and  Moin  to  provide  wall  stresses  in  a  trailing  edge  geometry.  It  is  given 
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by 

&  ,  ,  r\  d^i.RANS  o  (r  0r\ 

Q-V'  +  Vt) — —  =  °>  i  =  1»3)  l5-25) 

where  ul  is  the  RANS  eddy  viscosity  given  by  a  damped  mixing  length  model 

1±  =  l  +  Ky+  (l  —  e"^)2. 

The  boundary  conditions  are  taken  to  be 


m»,ransU  =  0) 

«t,RANs|m  =  ^i,LEs|m- 


Now,  the  terms  in  (5.17)  can  be  identified  as 


T>IZq(qr) 
Vq  q 


Vq 


q 


=  o, 
d 


/  .  r\^Ui,  RANS 


=  0. 


The  boundary  conditions  in  (5.18)-(5.19)  then  become 


h'q,qrq'r\m  ~  “i,RANsU> 
/r,bc(7>9r)  =  ^i,LEs|m- 


(5.26) 

(5.27) 


(5.28) 


(5.29) 

(5.30) 


The  adjoint  operator  can  be  determined  by  integration  by  parts  to  be 


du* 


i,RANS 


dy 


(5.31) 


Note  that  in  this  case,  vT  is  not  a  function  of  any  of  the  state  variables  so  it  can  be 
handled  directly  in  the  integration  by  parts.  The  boundary  terms  are 


^  f  f ,  ,  r\  (  *  ^u'h Rans  ,  Brans')  j 

BTr  -  j{v  +  vt )  ^u^rans — ^  wi,RANs  J  dxdy 


Vm 


(5.32) 


o 
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Given  these  equations,  it  is  now  possible  to  identify  the  RANS  sensitivities  at 
a  point  t/o  €  (0,  ym).  The  source  term  for  the  adjoint  equation  is  taken  to  be  the 
Dirac  delta  function,  8(x  —  x0),  in  both  the  u*  and  w*  equations  and  the  boundary 
conditions  are  chosen  to  be: 


ui,RANslo  —  ui,RANsU  ~  0- 


(5.33) 


With  these  values,  and  recalling  that  «(iRAnsIo  =  0,  (5.24)  becomes 


^rans(^o)  +  wRANS(x0)  —  /  / 

Jx  J  Y 


'“i.RANS 


d^j.RANS 


dy 


dxdz. 


(5.34) 


By  noting  that  the  RANS  equations  are  an  ODE  at  each  wall-normal  location,  the 
delta  function  used  as  the  RANS  adjoint  source  term  will  be  retained  in  the  x  and 
2  directions  in  the  adjoint  solution.  Hence  the  sensitivity  can  be  identified  as 


DUi, RANS^j^f) 
DiALEs(z,2/m,z) 


d^*,RANS 

dy 


(x,  Vm , 


*)> 


(5.35) 


Since  the  RANS  sensitivities  are  related  to  the  LES  state  at  only  one  point,  a  delta 
function  must  be  included  in  (5.7)  in  the  definition  of  the  operator  defining  the 
RANS  gradient  with  respect  to  the  LES.  Hence,  a  delta  function  must  similarly  be 
included  in  the  LES  adjoint  source  term  given  by  (5.9). 


5.3  Decoupling  the  Mean  Wall  Stress  from  the 
Control 

Several  different  control  formulations  proved  unsuccessful,  including  directly  replac¬ 
ing  the  fixed  target  in  Nicoud  et  al.  (2001)  with  the  RANS  velocity.  From  these 
results,  some  important  requirements  on  the  control  were  determined.  First,  the 
RANS  profile  must  be  able  to  correct  itself,  i.e.  if  there  is  an  error,  this  error  must 
decrease  with  time.  Second,  cost  functions  that  measure  a  ratio  of  the  velocity 
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profiles  tend  to  be  unstable  and  should  be  avoided.  Third,  in  order  for  the  control 
to  be  effective,  it  must  exert  most  of  its  effort  based  on  the  LES  state.  These  three 
results  imply  that  the  structure  of  the  cost  function  of  Nicoud  et  al.  (2001)  is  a 
good  one  to  use,  and  that  the  RANS  equations  must  somehow  be  coupled  to  mean 
flow  information,  in  this  case  either  the  mean  pressure  gradient,  or  equivalently, 
the  mean  wall  stress. 

Since  coupling  the  control  directly  to  the  RANS  equations  proved  ineffective, 
it  is  therefore  proposed  that  the  RANS  system  be  used  to  prescribe  the  mean 
wall  stress,  and  the  control  used  to  compute  the  fluctuations  about  this  mean. 
The  RANS  system  under  consideration  is  readily  adapted  to  that  purpose  since  the 
average  wall  stress  is  easily  obtained  by  using  the  system  as  an  algebraic  wall  model 
as  in  Wang  and  Moin  (2002).  The  mean  value  of  the  control  is  then  removed,  while 
all  the  fluctuations  are  retained.  This  can  be  formally  justified  by  considering  the 
control  not  to  be  the  local  wall  stress  but  instead  to  be  a  Fourier  coefficient  in  a 
Fourier  series  expansion  of  the  wall  stress,  i.e. 

Nx/2  Nz/2 

tw(x,  z)  =  J2  E  M"'1,  (5-36) 

n=-Nx/2  m=—Nz/2 

with  Nx  and  Nz  being  the  number  of  Fourier  modes  in  the  streamwise  and  spanwise 
directions,  respectively. 

In  this  formulation,  0nj7n  will  be  given  by  the  optimal  control  as  the  appropri¬ 
ate  Fourier  coefficients  for  all  n,m  except  for  n  =  m  =  0.  The  latter  coefficient, 
representing  the  mean  wall  stress,  will  be  prescribed  by  the  average  of  that  pre¬ 
dicted  by  the  algebraic  RANS-based  wall  model.  In  order  to  see  how  this  term  can 
be  removed  from  the  control  set,  consider  the  cost  function  gradient  identification 
equation  (4.24)  written  with  (5.36)  substituted  in  for 
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i(nNxx+mNzz) 


dt  r  r  /  Nl/ 2  Nz/ 2  - 

E  ^ 

Y  J  J  V  i=—Nx/2  j=-Nz/2 
Nx/2  Nz/2 

+w*  5Z  S  C,met(nNlI+mWz2)) 

i=-Nx/2j=-Nz/2 


which  can  be  rewritten  as 

W*/2  N*/2 

E  E 

n=-Nx/  2  m=-Nz/2 


DJ_~ 

Dcj) 


}=  E  E  fe»/  [u'eW-^^ixdz 

_ /vr  /o _ _ _  jit  /j  '  ^  ^ 

+<i>n,mj  J  w*ei{nNxX+mNzz)dxdz). 


By  applying  the  definition  of  the  Fourier  transform,  we  can  see  that 

Nx/2  Nz/2 

E  E 

n=-Nx/2m=-Nz/2 


(5.37) 


(5.38) 


(5.39) 


Thus,  the  gradients  can  be  identified  directly  from  the  Fourier  representation  of  the 
adjoint  solution: 


DJ 

DJ>1„ 

DJ 

D<t>ln 


=  U 


=  w 


n,— m* 


Since  the  gradient  of  each  Fourier  coefficient  can  be  identified  independently,  it  is 
possible  to  retain  only  a  subset  of  them  in  the  control  formulation.  Thus,  the  mean 
wall  stress  can  be  taken  out  of  the  control  set  and  the  optimal  control  problem 
can  still  be  solved.  It  should  be  noted  that  while  the  above  analysis  provides  a 
formal  argument  for  channel  flow,  a  more  general  Gram-Schmidt  procedure  could 
be  used  to  remove  the  mean  wall  stress  in  general.  All  that  is  required  that  it  can 
be  decomposed  into  mean  and  fluctuating  components.  In  more  complex  cases,  the 
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RANS  wall  stress  distribution  can  be  used  as  the  mean. 

Using  this  control  algorithm,  it  is  necessary  that  the  mean  wall  stress  is  deter¬ 
mined  by  RANS  at  the  old  time  to  prevent  the  controller  from  manipulating  the 
mean  stress,  while  the  fluctuating  component  can  be  determined  by  RANS  at  the 
new  time.  For  consistency,  however,  the  RANS  at  the  old  time  is  used  to  define 
target  in  the  cost  function. 

The  approach  can  then  be  summarized  as  follows:  first,  the  algebraic  wall  model 
of  Wang  and  Moin  (2002)  is  used  to  compute  the  mean  stress  at  the  wall.  In  this 
work,  the  velocity  inputs  into  this  model  come  from  the  velocity  at  the  first  LES 
grid  point  away  from  the  wall,  although  the  imposed  mean  pressure  gradient  ensures 
that  this  technique  using  the  input  from  any  point  would  correctly  predict  the  mean 
wall  stress  in  this  case.  In  order  to  compute  the  mean  wall  stress,  an  RANS  solution 
is  performed  in  the  wall-normal  direction.  These  RANS  velocities  are  then  used 
as  the  target  for  the  optimization  procedure.  By  combining  the  mean  wall  stress 
from  the  RANS  solution  and  the  fluctuating  wall  stress  from  the  control  routine, 
the  LES  wall  stress  is  found. 

Results  for  the  mean  velocity  profile  using  this  approach,  with  the  cost  function 
defined  over  the  first  three  grid  cells,  are  shown  in  Fig.  5.1  while  the  rms  profiles 
are  presented  in  Fig.  5.2.  Note  that  the  results  for  the  mean  velocity  profile  are 
almost  identical  to  what  was  obtained  by  Nicoud  et  al..  Thus,  the  model  retains  the 
fidelity  enabled  by  the  controller  but  is  predictive  in  that  no  a  priori  information 
was  needed  to  prescribe  the  target  profile. 

To  examine  the  robustness  of  this  controller,  it  is  tested  at  different  Reynolds 
numbers,  on  a  finer  grid,  and  with  different  SGS  models.  The  Reynolds  numbers 
considered  are  Rer  =  640,4000  and  20000.  First,  all  computations  are  performed 
on  a  grid  with  32  x  33  x  32  cells,  and  all  other  parameters  held  constant.  It  was 
found,  however,  that  the  convergence  rate  was  improved  at  ReT  =  640  by  increasing 
p.  This  is  likely  because  the  first  point  is  below  the  logarithmic  layer  and  so  more 
control  effort  is  required  to  increase  the  slope  of  the  mean  velocity.  For  all  Reynolds 
numbers,  though,  the  mean  velocity  profiles,  shown  in  Fig.  5.3,  accurately  capture 
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Figure  5.1:  Mean  velocity  profiles  at  ReT  =  4000, - :  u+  —  2.411og(y+)  +  5.2, 

- :  Piomelli  et  al.  (1989), - :  present  model, .  :  Nicoud  et  al.  (2001). 


the  logarithmic  profile. 

In  order  to  validate  that  the  method  is  independent  of  the  grid  spacing,  a  further 
computation  has  been  performed  at  ReT  =  4000  on  a  grid  with  64  x  65  x  64  cells. 
This  is  twice  the  resolution  in  each  direction  as  was  used  in  the  original  case.  All 
other  parameters  have  been  kept  constant.  The  mean  velocity  profile,  presented  in 
Fig.  5.4,  again  compares  favorably  to  the  logarithmic  profile,  demonstrating  that 
this  method  is  robust  with  respect  to  grid  and  Reynolds  number  changes. 

Another  important  robustness  issue  to  examine  is  the  sensitivity  of  the  control 
to  the  SGS  model.  Its  relevance  was  demonstrated  by  Baggett  et  al.  (2000)  in 
considering  the  LSE  wall  model  applied  in  different  numerical  environments.  While 
it  was  demonstrated  that  the  performance  of  the  LSE  approach  was  sensitive  to 
changes  in  the  numerical  method  and  grid  stretching,  the  technique  was  significantly 
more  sensitive  to  different  SGS  models  as  evidenced  by  the  large  over-prediction 
of  the  logarithmic  profile  when  the  Cabot  and  Moin  (2000)  procedure  was  used  to 
compute  the  Smagorinsky  coefficient  in  the  near- wall  region. 
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Figure  5.2:  RMS  of  velocity  fluctuations  at  ReT  —  4000, - :  Kravchenko  et  al. 

(1996), - :  Piomelli  et  al.  (1989), - :  present  model. 


In  this  study,  we  consider  three  models:  the  dynamic  Smagorinsky  model,  which 
was  used  to  obtain  the  previous  results,  the  mixed  similarity  model  Bardina  et  al. 
(1980),  and  the  procedure  of  Cabot  and  Moin  (2000).  The  stress  for  the  mixed 
similarity  model  is  given  by 

Tij  =  CSim(uiUj  —  UiUj )  +  DSM,  (5.40) 

where  the  model  constant  is  set  to  be  Csim  =  0.9,  (•)  denotes  the  same  test  filter 
as  used  in  the  dynamic  procedure,  and  DSM  stands  for  the  stress  coming  from  the 
dynamic  Smagorinsky  model.  The  results  using  these  three  models  in  conjunction 
with  an  algebraic  wall  model  Wang  and  Moin  (2002)  are  shown  in  Fig.  5.5.  It  is 
clear  that  these  SGS  models  can  produce  substantially  different  behavior  in  the  LES. 
However,  as  is  also  presented  in  Fig.  5.5,  when  they  are  used  with  the  control-based 
model,  the  near-wall  profiles  collapse  onto  the  logarithmic  profile  (again  without 
changing  the  optimization  parameters).  This  demonstrates  that  the  control  can 
account  for  different  SGS  modeling  errors,  while  the  LSE  wall  model  of  Nicoud  et 
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Figure  5.3:  Mean  velocity  profiles  computed  on  a  32  x  33  x  32  grid, - :  u+  = 

2.41  log(y+)  +  5.2, - :  ReT  =  20000, - :  ReT  =  4000, .  :  Rer  =  640. 


al.  Nicoud  et  al.  (2001),  based  on  one  instance  of  the  control,  could  not.  Further, 
note  that  away  from  the  wall,  the  quality  of  the  mean  profile  prediction  depends 
weakly  on  the  SGS  model  used.  It  can  be  seen  that  the  more  accurate  a  prediction 
the  model  makes  without  the  control,  the  more  accurately  the  mean  velocity  will 
be  captured  away  from  the  wall.  The  effect  of  the  control  is  to  ensure  that  the 
near-wall  mean  velocity  is  computed  correctly. 


5.4  Examination  of  the  Control  Efforts 

5.4.1  Correlations  Between  the  Control  and  Turbulent  Quan¬ 
tities 

Previous  efforts  to  understand  the  control  have  focused  on  its  effects  on  the  stress 
balance  and  turbulent  kinetic  energy  budget  (Nicoud  et  al.,  2001).  Such  an  ap¬ 
proach  aims  to  evaluate  the  control  based  on  the  changes  it  makes  to  the  time 


5.4.  EXAMINATION  OF  THE  CONTROL  EFFORTS 


85 


Figure  5.4:  Mean  velocity  profiles  for  ReT  =  4000, - :  vA  =  2.41  log(y+)  +  5.2, 

- :  64  x  65  x  64  cells, - :  32  x  33  x  32  cells. 


averaged  flow  field.  Here,  we  attempt  to  gain  insight  into  the  control  through  a 
different  statistical  measure:  its  correlations  with  turbulent  quantities.  These  cor¬ 
relations  allow  us  to  examine  what  the  control  is  reacting  to  in  the  turbulent  field. 
The  goal  is  to  better  understand  the  instantaneous  actions  of  the  control  rather 
than  the  changes  to  the  time-averaged  quantities.  In  addition,  these  results  may 
yield  insight  that  could  lead  to  feedback  models,  since  a  correlation  coefficient  with 
a  value  of  ±1  would  imply  that  a  perfect  feedback  controller  exists. 


In  this  investigation,  the  control  at  a  point,  4>{x,y ),  is  correlated  with  various 
turbulent  quantities  throughout  the  channel.  The  results  presented  this  section 
are  computed  at  Rer  =  4000  on  a  grid  with  32  x  33  x  32  cells  using  the  dynamic 
Smagorinsky  model.  If  h  is  some  function  of  the  turbulent  field,  h(x )  =  h(u,P ), 
we  can  compute  the  spatial  correlation  coefficient  with  the  expression 


((h(x)  -  S')  -  (f>)) 


(5.41) 
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Figure  5.5:  Effects  of  SGS  model  using  an  algebraic  wall  model  Wang  and  Moin 

(2002)  (left)  and  the  present  control-based  wall  model  (right),  -  :  u+  = 

2.41  log(y+)  +  5.2, - :  dynamic  Smagorinsky  model, - :  Cabot  and  Moin 

procedure  Cabot  and  Moin  (2000), .  :  mixed  similarity  model  Bardina  et  al. 

(1980). 


where  and  are  the  square  roots  of  the  variances  of  h  and  <j)  with  x  being 
a  location  in  the  channel  interior  while  x'  is  located  on  the  wall.  Averaging  is 
performed  over  wall-parallel  planes  and  in  time.  Because  of  the  spatial  homogeneity 
in  the  wall-parallel  plane,  the  correlation  is  only  a  function  of  the  spatial  separation 
C%(x  —  x').  The  results  for  the  maximum  correlation  coefficients  of  a  sample  of 
quantities  are  presented  in  Table  5.1.  In  all  cases  the  maximum  correlation  occurred 
in  the  second  wall-parallel  plane.  The  streamwise  and  spanwise  locations  of  the 
maxima  can  be  seen  in  the  figures  that  will  follow. 

These  results  indicate  what  the  control  is  and  is  not  reacting  to.  Perhaps  the 
most  interesting  result  is  the  lack  of  correlation  between  the  control  and  quantities 
related  to  the  shear  stress  balance  and  the  turbulent  kinetic  energy  (TKE).  The 
previous  work  of  Nicoud  et  al.  (2001)  focused  on  the  change  made  by  the  control  in 
decreasing  the  Reynolds  stress  and  increasing  both  TKE  production  and  dissipation 
near  the  wall  in  the  average  sense.  However,  these  results  suggest  that  the  control 
does  not  directly  respond  to  these  quantities,  as  their  correlation  coefficients  are 
quite  small.  Instead,  the  flow  is  manipulated  in  such  a  way  that  indirectly  changes 
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h 

max  |Q  | 

max  |Q | 

u 

0.46 

0.36 

V 

0.28 

0.19 

w 

0.19 

0.34 

TKE 

0.07 

0.03 

pk 

0.06 

0.03 

u'v' 

0.06 

0.03 

du/dy 

0.32 

0.22 

ux 

0.09 

0.29 

UJy 

0.17 

0.59 

Uz 

0.33 

0.21 

Ma 

0.09 

0.06 

du/dx 

0.59 

0.26 

Table  5.1:  Maximum  correlation  coefficients  for  the  streamwise  ( u )  and  spanwise 
(w)  boundary  conditions. 


these  flow  characteristics. 

To  examine  the  control  effects  in  more  detail,  it  is  useful  to  understand  their 
spatial  distribution.  Figures  5.6  and  5.7  show  both  the  horizontal  and  vertical  dis¬ 
tributions  of  the  coefficient  in  the  plane  of  its  maximum  value.  The  correlation 
coefficient  is  maximum  in  the  second  plane  from  the  wall  and  decreases  rapidly, 
consistent  with  the  near-wall  cost  function  construction.  Similarly,  the  effect  of 
the  control  is  also  local  in  the  wall-parallel  directions.  The  rapid  convergence  of 
the  steepest  descent  algorithm  is  likely  due  to  this  locality  since  the  optimization 
becomes  a  solution  of  many  local  problems  as  opposed  to  one  large  global  problem. 
Figures  5.6  and  5.7  also  indicate  that  the  wall  stresses  have  significant  structure. 
The  large  regions  of  high  positive  and  negative  correlations  aligned  with  the  stream- 
wise  axis,  in  the  case  of  (fiu,  and  in  the  spanwise  direction,  in  the  case  of  (f>w,  also 
suggest  that  velocity  gradients  will  be  more  highly  correlated  with  the  control. 

The  complexity  of  the  controller’s  actions  can  be  further  observed  by  considering 
the  spatial  distributions  of  additional  quantities  in  Figs.  5.8  and  5.9.  Patterns 
ranging  from  quadrupoles  to  “butterflies”  are  observed,  indicating  the  control  does 
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Figure  5.6:  Wall-parallel  spatial  distribution  of  the  correlation  coefficients  for  u !: 
coefficients  for  (f>u  (left)  and  <pw  (right)  at  y/h  =  0.09  (second  wall-normal  cell). 
The  scale  range  is  ±0.46  for  <pu  and  ±0.36  for  <pw,  with  white  being  a  large  positive 
value  and  black  a  large  negative  value.  The  domain  is  ±5  points  in  the  streamwise 
direction  and  ±6  points  in  the  spanwise  direction. 


Figure  5.7:  Wall-normal  spatial  distribution  of  the  correlation  coefficients  for  u 
coefficients  for  cf>u  (left)  and  <j>w  (right)  at  Az  =  0.  The  scale  range  is  ±.046  for  4>u 
and  ±0.36  for  <f>w ,  with  white  being  a  large  positive  value  and  black  a  large  negative 
value.  The  domain  is  ±5  points  in  the  streamwise  direction  and  5  points  in  the 
wall-normal  direction. 
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Figure  5.8:  Wall-parallel  spatial  distribution  of  the  correlation  coefficients  for  v: 
coefficients  for  4>u  (left)  and  <f>w  (right)  at  y/h  =  0.12  (second  wall-normal  cell). 
The  scale  range  is  ±0.28  for  0“  and  ±0.19  for  <j>w,  with  white  being  a  large  positive 
value  and  black  a  large  negative  value.  The  domain  is  ±5  points  in  the  streamwise 
direction  and  ±6  points  in  the  spanwise  direction. 


not  simply  apply  a  stochastic  force  to  the  flow.  In  all  cases,  the  control  is  only 
correlated  with  the  LES  over  a  few  grid  points,  reinforcing  that  the  controller  is 
solving  the  optimization  problem  locally.  The  highest  correlations  in  Table  5.1  are 
realized  with  spatial  distributions  involving  only  one  strong  peak,  as  shown  for  <f>u 
and  du/dx  in  Fig.  5.9  and  (f>w  and  uy  in  Fig.  5.10.  This  suggests  that  a  feedback 
controller  based  on  these  quantities  may  be  effective,  although  further  work  remains 
to  develop  a  feedback  law  independent  of  a  particular  simulation. 

These  correlation  results  have  implications  in  terms  of  applying  the  methodol¬ 
ogy  developed  here  to  more  complex  flows.  First,  the  fact  that  the  correlations  all 
peak  near  the  wall  demonstrates  that  the  near-wall  construction  is  indeed  appro¬ 
priate.  Tests  have  been  performed  using  a  cost  function  that  includes  the  entire 
domain  and  the  structure  of  the  correlation  coefficients  is  almost  identical  to  those 
shown  here.  In  addition  to  their  wall-normal  extent,  the  streamwise  and  spanwise 
correlation  lengths  suggests  that  local  averaging  may  be  sufficient  in  flows  with  no 
homogeneous  directions.  These  lengths  (measuring  approximately  ±3  grid  points 
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Figure  5.9:  Wall-parallel  spatial  distribution  of  the  correlation  coefficients  for 
du/dx:  coefficients  for  <fru  (left)  and  <f>w  (right)  at  y/h  =  0.09  (second  wall-normal 
point).  The  scale  range  is  ±0.59  for  4>u  and  ±0.26  for  <f>w,  with  white  being  a  large 
positive  value  and  black  a  large  negative  value.  The  domain  is  ±5  points  in  the 
streamwise  direction  and  ±6  points  in  the  spanwise  direction. 


Figure  5.10:  Wall-parallel  spatial  distribution  of  the  correlation  coefficients  for  u>y :■ 
coefficients  for  <f>u  (left)  and  <pw  (right)  at  y/h  =  0.09  (second  grid  cell).  The  scale 
range  is  ±0.17  for  4>u  and  ±0.59  for  <fiw ,  with  white  being  a  large  positive  value  and 
black  a  large  negative  value.  The  domain  is  ±5  points  in  the  streamwise  direction 
and  ±6  points  in  the  spanwise  direction. 
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in  each  direction)  can  be  used  to  define  the  averaging  operator  needed  in  the  cost 
function  formulation  while  still  supplying  the  controller  with  the  information  it 
requires. 

5.4.2  Structure  of  the  Near- Wall  Flow 

While  the  correlations  presented  in  the  previous  section  yield  some  information 
about  how  the  control  responds  to  the  flow,  it  is  also  useful  to  compare  instan¬ 
taneous  flow  realizations  between  controlled  and  uncontrolled  cases.  The  most 
striking  observation  is  that  very  little  appears  to  be  qualitatively  different  between 
simulations  using  the  algebraic  model  of  Wang  and  Moin  (2002)  and  the  control- 
based  framework  presented  here.  This  qualitative  comparison  is  demonstrated  in 
the  streamwise  velocity  fluctuations  at  the  first  wall-normal  plane;  contour  plots 
are  presented  in  Fig.  5.11  and  show  that  the  near- wall  structures  have  similar 
spatial  scales  and  organization  in  both  simulations.  This  indicates  that  much  of 
the  dynamics  occurring  near  the  wall  in  wall  modeled  simulations  is  qualitatively 
independent  of  the  wall  boundary  conditions.  Figures  5.12  and  5.13  show  that 
regions  of  high  streamwise  velocity  are  located  near  regions  of  intense  streamwise 
vorticity  and  wall-normal  velocity.  This  is  in  contrast  to  results  using  stochastic 
forcing  models  (Mason  and  Thomson,  1992;  Piomelli  et  al,  2003)  in  which  the  flow 
appears  much  more  de-correlated.  In  the  present  approach,  the  control  does  not 
randomly  force  the  flow  but  manipulates  it  specifically  because,  through  the  adjoint 
equations,  it  is  based  on  the  dynamics  of  the  simulation. 

Many  of  the  features  presented  in  Figs.  5.11  -  5.13  appear  qualitatively  sim¬ 
ilar  to  those  observed  in  resolved  lower  Reynolds  computations  (Moin  and  Kim, 
1982).  However,  the  structures  in  wall-modeled  computations  are  far  too  large  and 
do  not  accurately  represent  the  physics  of  near-wall  turbulence.  This  is  primar¬ 
ily  due  to  the  grid  spacing  as  the  spanwise  width  of  each  grid  cell  is  over  twice 
the  minimal  channel  spacing  required  for  self  sustaining  turbulence  (Jimenez  and 
Moin,  1991),  and  therefore  the  small  near-wall  dynamics  that  contribute  to  wall 
turbulence  are  not  resolved.  Since  the  structures  cannot  realize  their  correct  size, 
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Figure  5.11:  Contours  of  the  streamwise  velocity  fluctuations  at  the  first  wall- 
parallel  plane  with  control  (top)  and  without  control  (bottom).  Contour  levels 
are  from  —9 uT  to  13ur  for  the  controlled  case  and  from  —8 uT  to  10ur  for  the 
uncontrolled  case  with  dashed  lines  representing  negative  values. 
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Figure  5.12:  Contours  of  the  streamwise  vorticity  at  the  first  wall-parallel  plane 
with  control.  Contour  levels  are  from  —134 uT/h  to  203 uT/h  with  dashed  lines 
representing  negative  values. 

they  take  on  spanwise  widths  of  one  or  two  grid  cells.  Differentiating  the  field, 
to  determine  quantities  such  as  vorticity,  can  similarly  not  recover  the  underlying 
physics  because  the  grid  is  too  coarse  to  support  accurate  derivatives.  Rather, 
these  quantities  should  be  used  to  understand  the  dynamics  of  coarse  simulations 
since  these  derivatives  are  used  in  the  governing  equations.  Therefore  the  resulting 
structures  are  most  affected  by  the  grid  resolution  rather  than  the  wall  model. 

There  is  one  important  way  in  which  the  streamwise  velocity  fluctuations  do 
depend  on  the  control.  Without  control,  the  fluctuations  are  nearly  symmetric 
about  zero  with  a  maximum  amplitude  of  approximately  ±10uT.  When  the  control 
is  activated,  the  lower  bound  remains  constant  while  the  upper  bound  increases 
almost  40%.  This  affect  can  be  visualized  by  plotting  both  contours  using  the  scale 
of  the  results  without  control.  As  shown  in  Fig.  5.14,  this  scaling  exaggerates  the 
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Figure  5.13:  Contours  of  the  streamwise  vorticity  and  wall-normal  velocity  corre¬ 
sponding  to  the  line  in  Fig.  5.12.  Contour  levels  are  from  —134 uT/h  to  203 uT/h  for 
the  streamwise  vorticity  and  from  —3 uT  to  4 uT  for  the  wall-normal  velocity  with 
dashed  lines  representing  negative  values. 
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high-speed  regions.  While  the  actual  spatial  scales  of  these  regions  do  not  vary,  the 
control  increases  their  intensity  and  results  in  a  velocity  distribution  that  is  less 
symmetric  than  the  uncontrolled  case.  Because  of  the  underlying  dynamics  of  the 
simulation,  this  also  results  in  an  increase  of  intensity  of  the  streamwise  vorticity 
and  wall-normal  velocity. 

To  examine  this  aspect  of  the  control  in  more  detail,  comparisons  of  the  wavenum¬ 
ber  energy  spectra  of  the  streamwise  velocity  for  the  control-based  wall  model  and 
that  of  Wang  and  Moin  (2002)  are  presented  in  Fig.  5.15.  In  both  figures,  the 
spectra  are  normalized  by  the  square  of  the  streamwise  velocity  fluctuations  to 
determine  how  the  energy  distribution  changes  with  wave  number.  Very  little  dif¬ 
ference  is  observed  between  the  two  simulations  in  the  streamwise  spectra.  The 
spanwise  spectra  are  very  flat  compared  to  a  resolved  computation,  which  is  caused 
by  the  poor  resolution  of  this  grid.  There  is  apparently  a  shift  of  energy  from  low 
to  high  wavenumbers  caused  by  the  control.  In  addition,  the  two-point  correlation 
functions  of  the  streamwise  velocity  in  both  wall-parallel  directions  are  shown  in 
Fig.  5.16.  In  both  the  controlled  and  uncontrolled  cases,  these  correlations  are 
similar,  again  suggesting  that  the  control  does  not  change  the  size  of  the  near-wall 
features. 

The  data  implies  that  the  relevant  small  scale  in  these  simulations  is  the  grid 
spacing  rather  than  one  determined  by  the  physics.  To  investigate  this  issue,  a 
controlled  simulation  was  performed  on  a  grid  consisting  of  64  x  65  x  64  uniformly 
distributed  cells;  double  the  resolution  of  the  previous  simulations  in  each  direction. 
Figure  5.17  shows  the  contours  of  the  streamwise  velocity  fluctuations  in  the  first 
wall-parallel  plane.  The  structures  are  qualitatively  similar  to  those  on  the  coarser 
grid,  except  they  appear  finer.  This  observation  is  confirmed  by  the  streamwise 
and  spanwise  two-point  correlation  functions  of  the  streamwise  velocity,  presented 
in  Fig.  5.18.  In  both  directions,  the  structure  size  is  reduced  by  half  when  the 
grid  spacing  is  doubled.  Wall-normal  velocity  contours,  presented  in  Fig.  5.19, 
also  demonstrate  the  same  trend.  Note  that,  consistently,  the  size  of  the  structures 
increases  with  the  wall-normal  distance. 
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Figure  5.14:  Contours  of  the  streamwise  velocity  fluctuations  at  the  first  wall- 
parallel  plane  with  control  (top)  and  without  control  (bottom).  Contour  levels  are 
from  -8 uT  to  10itr  on  both  plots  with  dashed  lines  representing  negative  values. 
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Figure  5.15:  Energy  spectra  of  the  streamwise  velocity  in  the  streamwise  (left)  and 

spanwise  (right)  directions  at  the  first  wall-parallel  plane:  - :  no  control, - : 

control. 


Figure  5.16:  Two-point  correlation  function  of  the  streamwise  velocity  in  the 
streamwise  (left)  and  spanwise  (right)  directions  at  the  first  wall-parallel  plane: 
- :  no  control, - :  control. 
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Figure  5.17:  Contours  of  the  streamwise  velocity  fluctuations  at  the  first  wall- 
parallel  plane  with  control  on  a  grid  with  32  x  33  x  32  cells  (top)  and  on  a  grid 
with  64  x  65  x  65  cells  (bottom).  Contour  levels  are  from  — 9uT  to  13 uT  in  both 
cases  with  dashed  lines  representing  negative  values. 
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Figure  5.18:  Two-point  correlation  function  of  the  streamwise  velocity  in  the 
streamwise  (left)  and  spanwise  (right)  directions  at  the  first  wall-parallel  plane: 
- :  64  x  65  x  64  cells, - :  32  x  33  x  32  cells. 
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Figure  5.19:  Contours  of  the  wall-normal  velocity  in  a  streamwise  plane  on  a  grid 
with  64  x  65  x  65  cells.  Contour  levels  are  from  -4 uT  to  4ur  with  dashed  lines 
representing  negative  values. 


Chapter  6 

Conclusions  and  Future  Work 


The  two  primary  goals  of  this  work  were  to  reduce  the  computational  expense  of  the 
optimal  control  wall  stress  model  and  to  develop  a  predictive  model  based  on  this 
framework.  Both  of  these  goals  have  been  accomplished  for  high  Reynolds  number 
plane  channel  flow.  The  method  has  been  shown  to  produce  accurate  results  over 
a  range  of  Reynolds  numbers  and  is  grid  independent.  This  represents  a  significant 
advance  in  control-based  LES  wall  modeling. 

This  work  has  confirmed  the  results  of  previous  efforts  indicating  that  wall 
models  must  compensate  for  numerical  and  SGS  modeling  errors  that  are  present 
when  LES  is  performed  on  coarse  grids.  Many  earlier  wall  models  had  focused 
solely  on  compensating  for  the  unresolved  physics,  and  these  have  been  shown  to 
be  inadequate  on  a  coarse  grid.  A  wall  model  that  has  proven  able  to  compensate 
for  all  three  types  of  errors  is  the  control-based  wall  modeling  studied  in  this  work. 
The  earlier  version  of  this  method  was  unattractive  because  it  was  computationally 
expensive  and  not  predictive  due  to  the  fact  that  the  target  velocity  profiles  had  to 
be  prescribed. 

With  the  control-based  wall  model  developed  in  this  work,  the  computational 
expense  has  been  reduced  from  0(10)  times  that  of  the  LES  part  of  the  calculation 
to  0(1).  Reducing  the  expense  has  the  additional  effect  of  making  the  method  eas¬ 
ier  to  implement.  The  most  important  aspect  is  to  avoid  solving  the  optimization 
method  throughout  the  entire  domain.  It  has  been  shown  that  most  of  this  effort 
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is  wasted,  and  the  prediction  of  the  mean  velocity  profile  remains  accurate  when 
solving  both  the  adjoint  and  LES  equations  only  near  the  wall  in  the  optimization 
procedure.  Further  approximations  to  the  implicit  solutions  for  wall-normal  diffu¬ 
sion  and  the  pressure  have  also  contributed  to  cost  reduction  without  affecting  the 
accuracy  of  the  solution. 

The  final  aspect  of  reducing  the  computational  expense  is  to  reconsider  the 
derivation  of  the  adjoint  equations.  Previous  work  had  assumed  a  form  for  the 
adjoint  equations  based  on  a  Crank-Nicolson  numerical  technique.  There  is  no  a 
priori  justification  for  such  a  choice.  Often,  the  issue  is  whether  the  adjoint  equa¬ 
tions  are  derived  directly  from  the  discrete  system  or  instead  descretized  after  being 
derived  from  the  continuous  system.  While  the  latter  is  always  more  accurate,  when 
an  efficient  numerical  technique  is  required,  the  additional  expense  and  complexity 
of  such  an  approach  is  often  not  worth  the  modest  improvement  in  cost  function 
reduction  rate.  Since  we  are  already  approximating  the  cost  function  solution,  and 
since  the  actual  quantity  of  interest  is  the  time-averaged  mean  velocity  profile,  it  is 
reasonable  to  choose  the  adjoint  system  with  the  computational  cost,  rather  than 
maximum  accuracy,  as  the  primary  consideration. 

The  choice  of  the  discretization  is  tied  directly  to  the  definition  of  the  cost 
function.  It  has  been  demonstrated  that  defining  this  function  based  on  the  terminal 
time  rather  than  as  a  time  integral,  results  in  a  system  that  is  more  robust  with 
respect  to  temporal  discretization.  Thus,  the  adjoint  equations  are  discretized 
using  an  explicit  Euler  scheme.  By  using  such  a  strategy,  the  method  has  been 
made  significantly  more  computationally  tractable. 

The  other  outstanding  issue  of  previous  control-based  wall  models  is  the  reliance 
on  targets  chosen  a  priori,  making  them  not  predictive.  It  has  been  demonstrated 
in  this  work  that  velocity  profiles  obtained  from  RANS  can  be  used  to  determine  the 
targets  dynamically  during  the  simulation.  However,  coupling  the  RANS  solution 
to  the  LES  presents  several  difficulties.  The  primary  issue  is  that  the  controller  can 
manipulate  both  the  LES  and  RANS  solutions  if  they  are  coupled  incorrectly. 

By  coupling  the  RANS  with  the  controller  to  set  the  mean  wall  stress,  the  LES 
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and  RANS  solutions  can  be  made  compatible.  This  prevents  the  controller  from 
manipulating  the  mean  wall  stress  to  reduce  the  cost  function  by  producing  an 
unphysical  LES  computation.  It  has  been  shown  that  such  a  coupling  results  in  a 
simulation  as  accurate  as  those  which  use  a  fixed  target  profile. 

Given  that  a  wall  model  has  been  developed  which  is  computationally  simple, 
efficient,  and  predictive,  many  directions  for  future  work  are  available.  One  is  to 
apply  this  method  to  more  complex  flows.  By  using  more  general  near-wall  RANS 
treatments  and  coupling  them  to  the  controller,  the  techniques  presented  in  this 
work  can  be  extended  to  flows  in  complex  geometries.  Because  of  the  thin  nature  of 
the  wall  layer,  the  curvature  effects  can  often  be  neglected,  and  the  near- wall  region 
can  often  be  approximated  as  flat  plate  flow.  However,  the  averaging  operation  in 
complex  geometries  is  less  clear  since  there  is  often  no  homogeneous  direction.  The 
correlation  coefficients  presented  in  Chapter  5  demonstrate  the  the  control  only 
responds  locally  to  the  flow.  Therefore,  local  averages  should  be  used  to  define  the 
cost  function. 

Two  important  areas  for  future  work  are  that  of  heat  transfer  and  compressible 
flows.  To  extend  this  method  to  such  situations,  a  RANS  model  for  heat  transfer 
must  similarly  be  coupled  to  the  system.  Such  coupling  must  ensure  that  the 
controller  cannot  artificially  manipulate  the  RANS  heat  transfer  model,  as  this 
type  of  coupling  has  proved  detrimental  in  the  LES  momentum  equations.  A  similar 
approach  must  be  taken  for  compressible  flows,  with  RANS  models  providing  the 
mean  distribution  of  both  momentum  and  thermodynamic  variables  near  the  wall. 
The  control-based  wall  model  can  then  be  used  to  ensure  that  the  LES  solution 
matches  the  RANS  solution  near  the  wall.  The  mean  values  for  all  variables  based 
on  the  RANS  model  will  have  to  be  used  to  match  the  corresponding  LES  variables 
at  the  wall  to  correctly  couple  the  two  simulations. 

In  addition  to  these  areas,  further  work  will  be  needed  for  complex  physics 
simulations  where  the  fluid  mechanics  represent  only  part  of  the  problem.  One 
example  where  wall  modeling  has  already  been  utilized  is  in  aero-acoustics.  While 
acoustic  propagation  far  from  a  body  can  be  computed  using  the  Lighthill  analogy 
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with  the  solution  of  only  the  momentum  equations,  a  wall  model  coupled  with  a 
model  for  near- wall  pressure  fluctuations  could  be  used  to  obtain  even  more  accurate 
information  about  the  acoustic  field.  Another  interesting  application  would  be  to 
use  a  control-based  wall  model  in  flows  used  for  aero-optics  investigations.  In 
these  flows,  the  proper  affect  of  the  turbulence  on  the  optical  propagation  must  be 
modeled  near  the  wall  since,  as  was  shown  in  this  study,  the  fluctuations  in  this 
region  must  be  enhanced  to  obtain  a  correct  prediction  of  the  outer  flow.  Once 
this  model  is  applied  to  the  optics,  the  outer  flow  propagation  of  the  optical  beam 
can  be  handled  using  wall  modeled  LES.  Such  a  method  would  allow  aero-optical 
computations  to  be  performed  at  Reynolds  numbers  of  engineering  interest. 
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